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ABSTRACT 


A  comprehensive  computer  program  is  developed  for  two 
methods  of  spectral  analysis. 

The  Fast  Fourier  Transform  is  usea  as  the  basis  for  the 
modified  periodogram  method.  This  method  is  found  to  be 
highly  efficient.  It  is  used  as  the  primary  technique  for 
evaluating  Alnerta  climatic  data; 

The  Maximum  Entropy  Method  of  spectral  analysis  is 
investigated  and  found  to  be  lacking  m  ease  of  use  and 
accuracy.  This  method  is  intended  to  extract  very  long 
periods  present  in  a  data  set.  Tne  unavailability  of 
unbiased  data  casts  doubt  on  the  value  of  this  method  for 
climatological  investigations. 

An  evaluation  is  made  of  the  climatic  data  available  on 
microfiche  from  the  Atmospheric  Environment  Service, 
Fisheries  and  Environment  Canada. 

Eight  Alberta  stations  are  chosen  for  the  study.  A 
detailed  history  of  the  changes  of  the  meteorological  site 
is  compiled  for  each  station. 

A  five-year  periodicity  is  found  to  be  statistically 
significant  in  the  fall  temperatures  or  ail  eight  stations. 
Also,  three-  and  four-year  periodicities  are  found  to  be 
significant  in  the  winter  and  spring  temperatures, 
respectively,  of  a  majority  of  the  eignt  stations. 
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A 

period 

result 


thorough  investigation  is  made  of  the  five-year 
hut  no  evidence  is  found  to  indicate  that  it  could 
from  causes  other  than  true  climatic  fluctuations. 
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CHAPTER  1 


INTRODUCTION 

1  • 1  Periodicity  or  Fa n t a s y ? 

The  search  for  periodic  fluctuations  in  the  Earth’s 
climate  has  been  pursued  for  many  years.  The  motivation  for 
such  research  has  centred  around  tne  hypothesis  that  the 
climate  of  our  planet  may  repeat  itself  at  more  or  less 
regular  intervals.  Should  such  a  pattern  exist,  its 
discovery  would  contribute  greatry  to  long-range 
predictions.  However,  we  may  be  deluding  ourselves  at  times 
by  accepting  coo  quicxly  results  which  appear  at  first  sight 
to  be  "statistically  significant". 

A  favoured  approach  has  been  to  try  to  relate  the 
fluctuations  of  the  climate  to  the  fairly  regular 
periodicity  of  the  Zurich  sunspot  number.  A  brief  summary  of 
many  such  studies  was  compiled  by  Lawrence  (1965) .  Most  of 
the  studies  reviewed  in  this  paper  involved  graphical 
methods.  Some  studies  used  a  visual  correlation  technique 
which  is,  at  nest,  a  most  unscientific  approach,  while  other 
studies  considered  only  very  brief  periods  and  gave  no 
indication  of  the  extent  of  correlation  outside  such 
periods. 

Better  methods  of  analysis  have  been  developed 
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following  the  advent  of  modern  computers.  However,  this  has 
not  resolved  the  argument  over  the  sunspot  cycle  and  its 
relation  to  weather  and  climate.  A  good  example  of  this 
argument  is  shown  in  articles  written  by  Shapiro  (1975)  and 
Bain  (1976).  noth  authors  analyzed  a  315-year  record  of 
central  England  temperatures  compiled  by  Hanley  (1974). 
Shapiro  found  no  significant  evidence  of  sunspot-related 
periods  while  Bain,  using  a  different  analysis  technique 
called  the  Maximum  Entropy  Method,  found  evidence  of  the  22- 
23  year  double  sunspot  cycle.  The  method  used  by  Bain  is 
also  used  in  this  study.  However,  its  results  are  very 
difficult  to  assess,  as  will  be  seen  later,  and  should  not 
be  taken  at  face  value. 

Many  studies  in  which  time-series  analysis  was  used 
with  considerable  care  have  been  carried  out  over  the  last 
twenty  years.  A  brief  list  of  such  studies  would  include 
Landsbsrg  et  al.  (1  959),  Julian  (1971),  Joseph  (1973a, 
1973b) ,  Trenberth  (1975),  Bunting  et  al .  (1975), 
Bradley  (1976),  Jones  and  Kearns  (1976),  Madden  (1977),  Dyer 
and  Tyson  (1977)  . 


A  recent  article 

by 

Georgiades 

(  1977) 

has  provided 

evidence  of  a  four- 

to 

five-year 

cycle 

in  the  annual 

temperature  of  selected  stations  in  the  central  Canadian 
Prairies.  This  result  will  be  discussed  at  a  later  stage. 


1  • 2  The  Study 

The  present  study  attempts  to  combine  sound  statistical 
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methods,  recently  developed  methods  of  spectral  analysis, 
the  renewed  interest  in  climatic  change,  and  an  analysis  of 
the  climatic  conditions  over  one  of  the  more  important 
agricultural  regions  of  Canada. 

Based  on  the  numerous  studies  mentioned  in  the  previous 
section,  there  seems  to  be  little  hope  of  finding  evidence 
of  true  periodic  phenomena  in  temperature  and  precipitation 
data.  However,  a  study  of  the  distribution  of  variance  in 
the  frequency  spectrum  of  these  data  may  show  groups  of 
frequencies  wnich  contribute  significantly  to  the  total 
variance.  Conversely,  a  uniform  distribution  of  variance 
across  the  frequency  spectrum  would  suggest  that  there  is 
little  hope  of  using  climatological  information  for 
producing  long-range  forecasts. 

Care  must  be  taken  should  an  apparently  significant 
result  be  found.  Such  a  result  may  be  artificially  induced 
in  the  data  by  factors  such  as  the  data  processing  method, 
changes  of  location  of  the  observing  site,  aliasing  of 
frequencies  beyond  the  Nyquist  frequency,  or  leakage  of 
power  from  the  low  frequencies. 

To  date,  most  of  the  published  research  has  made  use  of 
the  well-known  autocovariance  method  of  spectral  analysis. 
However,  important  new  spectral  methods  have  been  revived  or 
developed  within  the  last  ten  years.  These  are  the  modified 
periodogram  analysis  using  the  Fast  Fourier  Transform  and 
the  Maximum  Entropy  Method  of  spectrar  analysis.  Each  of 
these  two  methods  has  distinct  advantages  over  the 
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autocovariance  method. 


The  Fast  Fourier  Transform  approach  produces  a 
considerable  saving  of  computational  time  while  providing  as 
much  detail  in  the  frequency  distribution  as  is  available 
when  using  the  ACY  method.  The  Maximum  Entropy  Method  offers 
the  potential  of  resolving  very  low  frequencies,  given  data 
which  contain  only  one  or  slightly  less  than  one  cycle  of 
the  pertinent  low  frequency.  Both  methods  will  be  discussed 
in  greater  detail  in  later  chapters. 

Rihiishi  (1976)  made  a  detailed  comparison  of  the 
modified  periodogram  and  autocovariance  methods  and  found 
that  the  two  methods  are  theoretically  equivalent.  However, 
the  modified  periodogram  permits  direct  calculation  of  the 
phase  angle,  while  the  autocovariance  method  requires  a 
separate  Fourier  analysis  in  order  to  derive  the  phase. 
Also,  his  calculations  of  the  computational  time  required 
were  based  on  a  comparison  between  a  heavily  smoothed 
periodogram  and  an  autocovariance  function  obtained  using  a 
fairly  small  maximum  lag.  Although  smoothing  is  applied  to 
the  periodogram  in  this  study,  the  unsmoothed  results  are  of 
equal  interest  since  they  can  be  analyzed  without  being 
concerned  about  the  effects  of  the  response  function  of  the 
filter.  A  significant  increase  of  computational  time  is 
required  to  Obtain  as  much  detail  from  the  autocovariance 
method  as  is  available  from  the  unsraoothed  periodogram. 
Thus,  the  modified  periodogram  was  preferred. 
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1.3  Weather  Parameters 


As  in 

any  study  of 

this  type. 

the  availability 

and 

homogeneity 

of  the  basic 

data  provide 

many  concerns. 

The 

basic  data 

collected  for 

the  study  are 

the 

published  monthly  mean  temperatures  and  the  total  monthly 
precipitation  amounts  (converting  snow  amounts  to  water 
equivalent)  at  eight  Alberta  stations.  The  selection 
procedure  is  described  in  Chapter  2  along  with  an  assessment 
of  accuracy  and  completeness . 

One  of.  the  major  problems  in  climatic  data  analysis  is 
inhomogeneity  resulting  from  changes  or  the  instrument  site. 
Some  of  these  changes,  such  as  installing  a  different  type 
of  thermometer  or  rain  gauge,  or  a  cnange  of  the  observer, 
are  generally  impossible  to  determine  prior  to  the  1950 ‘s 
since  comprehensive  records  of  such  changes  were  ill-kept 
before  that  time. 

Changes  of  location  or  of  the  surrounding  obstacles  are 
believed  to  be  of  greater  significance  and  can  usually  be 
ascertained  after  a  certain  amount  or  probing.  The  changes 
of  location  of  the  eight  stations  used  in  this  study  have 
been  determined  and  are  discussed  in  Cnapter  3. 

Different  groupings  of  the  data  are  considered  in  this 
study.  They  consist  of: 

1-  Using  all  values  at  one-montn  intervals  for 
each  parameter  to  determine  the  frequency 
distribution  for  frequencies  greater  than  one 
cycle  per  two  months.  This  analysis  requires  an 
appropriate  filter  in  order  to  suppress  the 
annual  cycle. 
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2-  Using  the  data  for  each  month  of  the  calendar 
year  in  order  to  determine  variations  in  the 
individual  months,  for  example,  to  determine 
the  frequency  spectrum  or  January  mean 
temperatures  or  January  precipitation  amounts. 

3-  Using  seasonal  means  in  order  to  determine 
variations  for  each  season.  The  frequent 
occurrence  of  convective  precipitation  in  the 
summer  in  Alberta  will  likely  make  this  season 
difficult  to  assess. 


A  comparison  of  results  for  each  station  should  also  be 
considered.  Such  an  analysis  may  shed  some  light  on  the 
areal  extent  of  significant  results,  i.e.,  it  may  indicate 
whether  apparent  cycles  are  strictly  a  local  phenomenon  or 
whether  they  can  be  identified  over  certain  areas  of  the 
province  if  not  over  all  stations  involved.  However, 
erroneous  conclusions  resulting  from  data  processing 


procedures  or  limited  sample  sizes  common  to  all 


would  not  necessarily  be 


revealed  by  such  analyses. 
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CHAPTER  2 


DATA 


2 .  1  Selection 

Within  the  province  of  Alberta,  there  are  a  number  of 
locations  which  have  weather  records  covering  periods  of  the 
order  of  70  to  100  years. 

The  first  step  was  to  determine  which  weather  or 
climatological  stations  had  kept  records  that  combined  the 
following  criteria: 

1-  Length,  of  records  exceeding  50  years. 

2-  Complete  temperature  and  precipitation  records. 

3-  Insignificant  or  no  change  of  .location  of  the 
observing  site  during  the  length  of  the  record. 

The  most  convenient  source  that  provided  this 
information  is  the  Climatological  Station  Data 
Catalogue  (1976).  A  total  of  27  stations  was  chosen  as  a 
possible  set  of  benchmark  stations  for  the  project;  very  few 
of  these  27  stations  met  all  three  criteria  but,  of  the  over 
600  stations  listed  in  the  catalogue,  these  27  stations  came 
closest  to  meeting  the  criteria.  As  will  be  explained  in 
Chapter  3,  this  catalogue  is  not  a  fully  reliable  source  of 
information  nut  it  does  provide  an  excellent  guide  for 
making  an  initial  choice. 
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The  next  step  was  to  examine  the  records  available. 
This  was  readily  done  using  the  microfiche  station  records 
available  at  the  Western  Regional  Office  (WRO)  of  the 
Atmospheric  Environment  Service  (AES)  located  in  Edmonton. 

The  WRO  has  available,  on  microfiche,  the  complete 
monthly  record  of  every  station  in  Alberta  up  to  December 
1972.  Each  microfiche  plate  provides  such  items  as  monthly 
mean  maximum,  mean  minimum  and  mean  temperature,  extreme 
temperatures  for  the  month,  total  monthly  precipitation, 
rainfall  and  snowfall,  and  many  other  monthly  values. 

Of  interest  were  the  monthly  mean  temperature  and  the 
total  monthly  precipitation.  A  scrutiny  of  the  microfiche 
records  provided  immediate  and  detailed  information  as  to 
the  completeness  of  the  records.  This  process  reduced  the 
total  number  of  suitable  stations  to  19. 

Copies  of  the  pertinent  information  were  obtained  for 
these  19  stations.  Values  for  the  period  after  1972  were 
obtained  from  the  Monthly  Record  (1973-1975),  published  by 
the  AES. 

Of  the  19  stations,  8  were  chosen  as  the  main  stations 
to  be  used  in  the  project.  They  are: 

1-  BEA  -  Beaverlodge  CD  A  (Canadian  Department  of 
Agriculture)  -  60  years  (of  records) . 

2-  YXD  ~  Edmonton  City  -  93  years. 

3-  FVR  -  Fort  Vermilion  CDA  -  68  years. 

4-  LAC  -  Lacombe  CDA  -  68  years. 

5-  LTH  -  Lethbridge  CDA  -  68  years. 


6-  WAN  -  Manyber  r.ies  CD  A  -  4  7  years. 

7-  YXH  -  Medicine  Hat  -91  years. 

8-  RAN  -  Ranfuriy  -71  years. 

Although  MAN  falls  short  of  the  50-year  criterion,,  it 
was  chosen  Decause  of  the  completeness  of  its  record;  the 
value  of  such  a  record  will  become  evident  later  on. 


2.2  Verification 

In  order  to  determine  the  accuracy  of  the  data  kept  on 
microfiche,  a  brief  comparison  of  these  data  was  made  with 
the  published  values  in  the  Monthly  Record  (1880-1972) .  This 
comparison  revealed  a  number  of  discrepancies.  As  no  one  at 
the  WHO  could  provide  definite  information  concerning  the 
accuracy  of  the  microfiche  data,  a  comprehensive  comparison 
of  data  was  undertaken  to  determine  the  extent  of  these 
discrepancies.  Where  available,  daily  values  were  used  to 
recalculate  the  monthly  values  when  these  values  did  not 
agree  with  the  microfiche  record. 

After  a  number  of  discrepancies  were  found,  a  letter 
was  sent  to  the  Information  Section  of  the  Headquarters  of 
the  AES  in  Toronto,  Ontario,  requesting  verification  of  the 
microfiche  data  in  cases  where  it  did  not  agree  with  the 
published  monthly  values.  The  comparisons  that  were 
undertaken  and  the  results  are  shown  in  Table  2.1.  It  should 
be  noted  that  tolerances  of  .2F  («11C)  for  temperature  and 
.02  inches  (.5mm)  for  precipitation  were  allowed  since 


earlier  data  would  have  been  recalculated  by  computer 
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resulting  in  slight  differences  due  to  round-off  of  daily 
and  monthly  values. 

Table  2.1  shows  that  of  350  discrepancies  found  by 
comparing  8565  values,  88  were  verified  and  all  findings 
supported  the  values  given  in  the  microfiche  record.  Of  the 
350  discrepancies,  207  were  for  differences  greater  than  one 
degree  Eahrenneit  (.56C)  or  one  tenth  of  an  inch  (2.5mm)  of 
precipitation;  61  of  the  88  values  resolved  were  also  in 
this  category,  with  some  differences  being  as  large  as  one 
order  of  magnitude  or  more.  Also,  103  of  the  350 
discrepancies  occurred  in  the  very  early  data  (prior  to 
1891).  The  breakdown  of  these  103  values  is  contained  in 
Table  2.1,  where  the  values  in  parentheses  represent  the 
period  prior  to  1891  and  are  included  in  the  number  given 
for  the  total  period  of  comparison. 

As  a  result  of  these  findings,  dr.  Derek  Aston  was 
contacted  in  Toronto.  Ke  explained  the  procedure  used  to 
compile  the  microfiche  records  which  are  available  for  most 
stations  across  Canada. 

The  process  began  in  the  early  1960‘s  and  required 
several  years  to  complete.  Trained  Meteorological 
Technicians  verified  the  original  weather  records  for  each 
station  up  to  the  time  when  comprehensive  verification  of 
data  became  routine.  Errors  and  conflicts  found  in  the 
original  records  were  resolved  using  all  available 
information  and  the  technicians'  knowledge;  errors  and 
conflicts  whrch  could  not  be  resolved  satisfactorily 
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resulted  in  a  missing  value  on  tie  microfiche.  Major 
differences  between  the  new  computed  values  and  the  original 
published  information  were  again  verified.  The  objective  of 
this  program  was  to  produce  data  whrch  had  been  processed 
uniformly  for  the  entire  length  of  the  record  and  the  final 
product  was  made  available  on  microficne. 

Although  there  exists  a  slight  possibility  that  errors 
may  be  present  in  the  microfiche  record  due  to  bey-punch 
errors,  this  author  feels  that  the  findings  of  the 
comparison  carried  out  reduce  the  frequency  of  such  errors 
to  insignificant  proportions.  The  time  and  effort  involved 
in  compiling  the  microfiche  record  have  resulted  in  one  of 
the  most  accurate  amd  homogeneous  sets  of  weather  records 
available. 

2 . 3  Summary 

Having  accepted  the  microfiche  records,  the  next  step 
was  to  determine  values  for  the  few  missing  entries. 

In  carrying  out  the  comparison,  some  of  the  values 
entered  as  missing  on  the  microfiche  were  found  in  the 
Monthly  Record.  Some  of  these  missing  values  were  also 
investigated  by  Mr.  Aston.  However,  no  changes  were  made  to 
the  microfiche  record,  indicating  that  the  available  data 
for  that  montn  did  not  meet  specified  requirements. 

First-guess  estimates  of  the  missing  values  were  made 
using  the  values  published  in  the  Monthly  Record  or,  if  also 
missing  in  this  source,  by  estimating  a  value  based  on  other 
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Table  2.2.  Summary  of  data  sets.  (T  -  Temperature*. 

P  -  Precipitation.) 


Data 

Set 

Period 

Total 

Number 

Number 

Number 

Number 

Miss ing 

Found* 

Est. 

BEA 

T 

Sep/ 191 5- Dec/ 1975 

724 

0 

0 

0 

P 

Apr/ 19  15-Dec /1975 

729 

2 

1 

1 

YXD 

T 

Sep/1881 -Dec/ 19 75 

1  132 

4 

0 

u 

P 

Apr/ 1 8 82- Dec/ 1975 

1  125 

1 

0 

1 

FVR 

T 

Jul/1 9  08-Nov/l 9  75 

809 

0 

0 

0 

P 

Jul/1 9  08 “Nov/ 19  75 

809 

5 

1 

4 

LAC 

T 

Dec/ 19  0  7-  Dec/ 1  9  75 

817 

1 

1 

0 

P 

Mar/1 9  08-Dec/l 9  75 

814 

6 

3 

3 

LTH 

T 

Mar/1908-Dec/1975 

814 

0 

0 

0 

P 

Feb/ 1 9  08 -Dec/ 1 9 7  5 

815 

0 

0 

0 

MAN 

T 

Sep/ 19 28 “Dec/1 9  75 

568 

0 

0 

0 

P 

May/ 19 28- Dec/ 19 75 

572 

2 

0 

2 

YXH 

T 

J  ul/ 1 8  8  4 - D  ec/ 19  75 

1  098 

4 

1 

3 

P 

Oct/ 1 8  83 - D  ec/ 19  7  5 

1  107 

15 

11 

4 

RAN 

T 

Jan/ 19  05-Dec/ 19 75 

852 

0 

0 

0 

P 

Jan/1905-Dec/1975 

852 

5 

4 

1 

Sub 

T 

6814 

9 

2 

7 

Total 

P 

6823 

36 

20 

16 

Total 

13637 

45 

22 

23 

*  Data 

•Found*  refers  to  data  published 

in  the  Monthly 

Record  but  rejected  when  the  microfiche  data  were  compiled. 

stations 

in  tne  same  area. 

Table  2.2 

snows 

the  total 

periods 

of  data 

involved  and  the  numbers  of 

missing 

data  values. 

Only  45  values  were  missing  out 

of  a 

total  of 

13,637 

from 

the  eignt  stations  chosen,  and 

34  of 

these  occurred  in 

1916 

or 

earlier.  Thus,  for 

t he  period  19  17 

to  1975, 

only  11 

values  were  massing  from  a  total  or  11,050  values,  and  seven 
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of  these  were  found  in  the  Monthly  Record. 

Further  discussion  of  the  missing  values  will  be  made 
in  Chapter  7.  It  suffices  to  say  at  this  time  that  the  data 
prepared  for  this  study  are  believed  to  be  as  sound  as 
possible. 
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CHAPTER  3 


STATION  HISTORIES 

3. 1  Sources  of  Informa tion 

Any  statistical  analysis  of  meteorological  information 
must  take  into  account  changes  of  location  of  the  observing 
sites.  The  listings  in  the  Cliraatoiogical  Station  Data 
Catalogue  provide  a  history  of  each  station  in  terms  of 
latitude,  longitude  and  elevation.  However,  comprehensive 
on-site  inspection  by  AES  of  observing  stations  did  not 
begin  until  tne  early  1950* s;  thus  earlier  knowledge  of 
station  locations  is  often  very  meager  and  subject  to  much 
doubt.  An  appeal  for  improvement  was  recently  made  by 
Catchpole  and  Ponce  (1976). 

Since  the  Catalogue  showed  no  change  in  the  location  of 
Medicine  Hat  since  1883,  doubt  was  immediately  cast  on  the 
accuracy  of  its  listings;  it  did  not  seem  likely  that  the 
present  Medicine  Hat  Airport  would  be  located  at  the 
original  observing  site. 

In  order  to  eliminate  as  much  of  this  doubt  as 
possible,  a  tnorough  search  of  all  information  available  was 
made.  The  nest  local  sources  of  information  are  the 
inspection  report  files  at  the  WRO.  These  files  provided 
good  histories  of  Medicine  Hat  and  Runfurly.  A  history  for 


15 


4 . 


_ 


*  *  «* 


•t  /  •  - 

V* *  »• 

V 

• 

t 

«  >  • 

16 


Edmonton  was  readily  available  in  a  separate  source.  As  for 
the  CDA  stationsr  little  information  could  be  gained  from 
the  files  so  that  letters  were  sent  to  these  places 
requesting  historical  information  if  available.  All  replies 
were  very  complete  and  most  included  maps  showing  the  exact 
changes  of  the  meteorological  sites. 

The  history  of  Medicine  Hat,  as  found  in  the  WHO  files, 
was  not  totally  clear.  There  were  references  to  commercial 
establishments,  some  of  which  are  no  longer  in  existence.  A 
visit  to  the  Alberta  Government  Telephones  archives  was 
necessary  to  determine  the  exact  locations  of  the 
meteorologicar  sites. 

The  remainder  of  this  chapter  outlines  the  history  of 
each  station.  The  possible  consequences  of  observing  site 
changes  will  be  discussed  along  with  the  numerical  results 
in  Chapter  7. 

Appendix  A  contains  a  topographical  map  for  each 
station  showrng  the  various  positions  of  the  observing 
sites. 


3.2  Beaverlodge  CDA 

A  letter  received  from  Mr-  H.  Hall,  Agricultural 
Meteorologicar  Technician  at  the  experimental  farm  provided 
detailed  information  on  the  Beaverlodge  site. 

The  first  site,  at  an  elevation  or  750ra  (2460  feet), 
was  to  the  North  of  the  present  site  at  an  elevation  of  739m 
(2425  feet);  the  distance  between  the  two  sites  is  373.7m 
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1 1 226  feet)  .  The  terrain  in  this  area  is  slightly  rolling. 

Mr.  Hall  also  advised  that  prior  no  the  official  move 
of  the  observing  site  in  1958,  a  three-year  comparative 
study  was  made.  Although  the  original  records  could  not  be 
found,  this  study,  referred  to  by  carder  (1962),  will  be 
discussed  in  Chapter  7. 

The  official  position  of  the  present  site  is  55°128 
North  119°25‘  West  at  an  established  elevation  of  731.5m 
(2400  feet)  anove  sea  level  (ASL)  .  (Established  position  and 
elevation  may  differ  from  actual  values  because  of  the  scale 
of  the  topographical  map  used;  WHO  uses  maps  with  a  scale  of 
1/250,000.) 


3 . 3  Edmonton  city 


The  hist or 
Annual  Meteoro 
to  be  used  star 


y  of  Edmonton  was  readily  available  in 
logical  Summary  for  Edmonton  (1974) .  As 
t  in  September  1881,  tne  history  prior 


the 

data 

to 


this  date  is  presented  for  completeness  only. 


The  following  text  is  taken  from  the  Annual  Summary: 


History 

On  Decembe 
telegrapher,  c 
meteorological 
Instruments  incl 
one  maxrmum  an 
gauge,  one  (comm 
and  a  .Forsters 
vane.  On  May  19, 
No. 2389,  was  fir 
Dairy  obse 
with  minor  inter 
repairing  the 
anemometer  was  s 
several  evening 


r  29,  1879,  Mr.  G.S.  Wood,  a 

ompleted  installation  of  a 
station  at  Fort  Edmonton, 
uded  two  standard  thermometers, 
d  one  minimum  thermometer,  a  rain 
on)  thermometer  screen  and  shed. 
Small  anemometer  with  double  tail 
1880,  a  Greens  Small  Barometer, 
st  used. 

rvations  commenced  July  11,  1880, 
ruptions  when  Mr.  Wood  was  absent 
telegraph  line.  Evidently,  the 
et  up  outside  the  Fort,  because 
observations  of  the  wind  are 


18 


missing  with  the  notation:  "Fort  gate  locked 
prematurely;  could  not  get  our  to  take  wind 
gauge" . 

In  1882,  the  meteorological  station  was 
turned  over  to  the  Hudson  Bay  company.  Hr.  C.F. 
Young  and  his  brother  took  the  reports  at  a  site 
just  north  of  Jasper  Avenue  on  what  is  now  115th 
Street.  Ihis  is  just  3.2km  (2  miles)  south  ox  the 
present  weather  station  and  is  about  the  same 
altitude. 

In  April  1912,  Mr.  S.M.  Holmden  became 
responsible  for  the  observations  and  the  site  was 
changed  to  63rd  Street  in  the  Higniands.  Altitude 
of  this  station  was  only  7.6m  (25  feet)  lower  than 
the  present  Airport  site. 

Mrs.  E.  Owen  took  over  the  observing  duties 
in  1915  and  continued  in  this  post  till  the 
climatological  station  was  closed  i.n  1942. 
(Meteorological  reports  for  Edmonton  use  the  data 
collected  at  the  Airport  since  1937.) 

Present  location 

The  approximate  location  of  the  Airport  is 
53°35*  North,  113°30®  West.  The  orficial  elevation 
of  the  Airport  is  670.5m  (2200  feet).  The  Airport 
is  surrounded  by  the  City  of  Edmonton  and  is 
approximately  a. 8km  (3  miles)  northwest  of  the 
city  centre.  The  meteorological  office  was  first 
moved  to  the  Airport  as  a  permanent  arrangement 
approximately  the  first  of  September,  1937.  it  was 
then  located  in  what  is  now  the  Western  Air motive 
Hangar.  Approximately  the  first  or  November,  1942, 
it  was  moved  to  the  Airport  Administration 
Building . 


A  new  terminal  building  has  ne 
Weather  Office  was  relocated  early  in  19 
base  ends  with  December  1975  and  is  not 


en  completed  and  the 
76;  however  the  data 
affected  by  this 


latest  move. 


3. 4  Fort  Vermilion  CPA 

A  letter  received  from  Mr.  B.  Siemens,  the 
Superintendent  of  the  experimental  farm,  relayed  the 
following  information: 
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From  1908  until  the  winter  of  1935-36,  the 
site  was  located  on  a  river  flat  on  the  south  bank 
of  the  Peace  River.. . .  at  an  estimated  elevation 
of  259m  (850  feet) .  I  beiieve  that  this 

meteorological  site  was  in  an  open  area  well  away 
from  busn. 

Since  1936,  the  raeteorologica 1  site  has  been 
on  the  Experimental  Farm,....  at  an  elevation  of 
277ra  (910  feet).  The  site  is  in  an  open  field. 

The  move  was  approximately  7.65km  (4.75  miles) 

eastwards.  The  present  location  of  the  station  is  58°23c 

North,  116°021  West,  and  the  established  elevation  is  279m 

(915  feet)  ASL. 

3 • 5  haccrabe  CD A 

A  letter  received  from  Mr.  J.R.  Gillespie,  the 
Meteorological  Observer  at  the  farm,  advised  that  the 
observing  site  had  been  moved  only  once  in  its  history. 

The  site  was  moved  1.15km  (0.72  mile)  south westwards  to 
a  more  open  area.  This  move  was  made  on  July  1,  1972. 

The  terrain  in  this  area  is  very  revel.  Thus  the  move 
resulted  in  a  negligible  change  of  elevation  (approximately 
one  meter  lower  than  the  previous  location) . 

The  official  position  of  this  site  is  52°28‘  North, 
113°45*  West,  and  its  established  elevation  is  8u7m  (2780 
feet)  ASL. 

3.6  Lethbridge  CD A 

A  letter  was  received  from  Mr.  E.H.  Hobbs,  a  research 
scientist  at  the  CDA,  outlining  the  history  of  the 


meteorological  site. 


20 


Since  observations  began  in  1908,  the  site  has  been 
moved  only  once  -  in  June  of  1966.  The  move  was 
approximately  602m  (1975  feet)  so uth west  wards,  with 
negligible  change  in  elevation  (a  decrease  of  3.4m  over  very 
flat  terrain) . 

The  present  location  is  49°42f  North,  112°47*  West,  and 
the  established  elevation  is  899m  (2950  feet)  ASL. 

3.7  da n_yberri.es  CPA 

flanyberr-A.es  CDA  is  a  substation  of  Lethbridge  CD  A. 
dr.  E.H.  Hobbs,  in  his  letter  about  Lethbridge  CDA,  outlined 
the  history  of  the  meteorological  station  at  Hanyberries 
CDA. 

Since  observations  began  in  1928,  the  site  has  been 
moved  twice,  both  moves  were  done  because  the  shelterbelt 
had  grown  too  high  over  the  years. 

The  original  location  was  a  few  nundred  meters  west  of 
the  Headquarters  Building  (HQB) .  In  1944,  the  site  was  moved 
to  a  location  129.5m  (425  feet)  east  of  the  HQB;  and  in 
1949,  the  site  was  moved  to  a  location  259m  (850  feet) 
northnortheas t  of  the  HQB  resulting  in  an  elevation  increase 
of  7.6m  (25  feet)  . 

The  present  location  is  49°07'  North,  110°28'  Rest,  and 
the  established  elevation  is  934m  (3065  feet)  ASL. 

3.8  Medicine  Hat 


A  newspaper  article  written  in  1933  by  an  unnamed 
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author  was  round  in  the 
concise  history,  taken  from 
1956,  was  also  found  in 


inspection  files  of  the  WHO.  A 
this  artrcle  and  updated  to 
the  files.  The  history  is  as 


follows: 


The  first  month  for  which  observations  are 
complete  is  that  of  September  1883.  The  observer 
was  hr.  J.L.  Ewart.  The  station  was  situated  on 
Second  Street  on  the  site  now  occupied  by  Taylor 
Bros,  (in  1933).  The  elevation  was  665.3m  (2183 
feet)  above  mean  sea  level. 

In  1891,  the  station  was  taken  over  by  Mr. 

J.K.  Drinnen,  and  moved  across  tne  street  to  a 
point  where  the  Cawker  Drug  Store  now  stands. 

Another  change  was  made  rn  1901,  when  Mr. 
Waiter  crosskill  took  over  as  observer.  The 
station  was  located  where  the  Moore  furniture 
store  is  now  situated  (about  90m  south  of  City 
Hail) . 

In  -January  1911,  Mr.  H.H.  Hassard  was  the 
observer.  The  station  was  moved  to  the  corner  of 
Fourth  Street  and  Ash  Avenue,  where  observations 
were  taken  until  Mr.  Hassard* s  death  in  1929.  Mr. 

W.H.  Phelps  became  acting  observer  until  August 
1930. 

In  August  1930,  Mr.  C.  Pickering  was 
appointee  full  time  weather  observer  to  set  up  and 
operate  the  weather  station  in  connection  with  the 
opening  of  the  Prairie  Air  Mail  Service.  The 
station  was  set  up  temporarily  at  the  High  School 
site  where  the  elevation  was  70b. 5m  (2318  feet), 
then  moved  to  the  airport  at  an  elevation  of 
720.8m  (2365  feet)  in  the  Spring  of  1931.  When  the 
Prairie  Air  Marl  was  abandoned  in  March  1932,  the 
weather  station  was  moved  back  to  the  city  to  a 
location  on  Crescent  Heights  ar  an  elevation  of 
703.  1m  (c307  feet)  . 

In  1933,  the  station  was  moved  to  a  location 
immediately  west  of  Division  Avenue  on  Third 
Street  at  an  elevation  of  704m  (2310  feet) . 

The  station  was  moved  back  to  the  Airport 
Radio  Range  station  in  1938  at  a  location  on  the 
north  side  of  the  Airport  which  is  located  about 
3.2km  (2  miles)  southwest  of  the  city  centre. 

There  are  no  records  of  any  changes  after  1938.  In 

brief,  the  observing  site  was  located  in  the  South 

Saskatchewan  River  valley  from  1883  to  1930,  and  it  has 
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remained  out  of  the  valley  since  1931. 

.  The  present  location  is  50°01*  North,  11C°43'  West,  at 
an  established  elevation  of  720. 8m  (23o5  feet}  AS1. 

3 . 9  Ran fur ly 

A  brief  history  found  in  the  inspection  files  at  the 
WHO  shows  that  the  climatological  program  has  been 
maintained  by  the  same  family  since  the  first  of  January, 
1905.  Presently,  the  grandsons  of  the  original  observer  are 
performing  the  duties. 

The  instrument  site  was  moved  about  C.8km  (.5  mile) 
westwards  in  the  Fall  of  1931.  Previous  to  this,  the 
location  was  closer  to  brush  and  more  sheltered.  From  the 
fall  of  1931  to  June  18,  1951,  the  thermometer  screen  was 
situated  in  a  more  sheltered  location,  some  12m  (40  feet)  to 

the  east  of  its  present  site.  The  move  in  1931  resulted  in 
no  appreciable  change  in  elevation. 

The  position  of  the  station  is  5o°27*  North,  111°399 
West,  and  the  established  elevation  is  685.8m  (2250  feet) 


ASL 


CHAPTER  4 


EAST  FOURIER  TRANSFORM 


4 . 1  The or  y 

The  Fast  Fourier  Transform  (FFT) ,  first  developed  in 
the  early  1900‘s,  returned  to  prominence  in  the  mid  1960 8 s 
as  a  method  capable  of  saving  considerable  computational 
time  (Cochran  et  al.  ,  1967;  Bingham  ec  ai„ ,  1967;  Cooley  et 
al.  ,  1967)  . 

The  basic  formula  for  the  transform  is 


where  xK  represents  a  time  function,  y-  represents  the 

transform,  ana  j=0 ,  1, . . . . , g-  1 . 

Equation  (4.1.1)  forms  the  basis  of  the  FFT  even  though 
it  is  the  definition  of  the  Discrete  Fourier  Transform 

(DFT)  . 

In  calculations  of  the  FFT,  the  data  set  is  factored 
into  multiple  components  of  two  before  using  (4.1.1).  Using 
matrix  notation,  the  procedure,  as  described  by 
Robinson  (1967),  consists  of  taking  two  Ixq  row  vectors  X 
and  Y  and  a  qxq  matrix  K  to  obtain 


Y  =  XW 


(4.  1.2) 


where 
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*  =  t  y» » y, 


W  =  [e 


-iZITk 


JA] 


k, j  =  0 ,  1 , 2,  . . .  .  ,  g-  1 


The  matr-LX  W  is  symmetric  and  has  the  property  that 

WW*r  =  ql.  (4 .  1.3) 

The  FFT  algorithm  factorizes  the  matrix  W  into  n  sparse 


matrices  S, , S2 ,  . . • . , S*  and  a  permutation  matrix  P  such  that 


W  =  S,  S2 . Sn  P  (4.  1. 4) 

where  n  =  log^g. 

The  sparse  and  permutation  matrices  can  be  illustrated 
for  q=23=8.  If  we  define  w  as 

«  = 


then  the  matrix 

}J  can 

be 

written 

as 

VJ 

=  [  w 

J 

k,  j= 

0,1,2 

«.  *0  O  <  ^  p  ^ 

Thus  we  nave 
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The  order 

of 

the 

powers 

of 

W 

in 

the 

sparse  and 

permutation 

matrices 
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the 
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This  yields 


P  = 


c  = 

I 


S2  ~ 


and 
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Thus  the  algorithm  requires  tnat  the  number  of  data 
points  be  a  power  of  two.  If  the  number  of  values  in  the 
data  set  is  not  a  power  of  two,  zeroes  must  be  added  until 
the  total  number  q  meets  this  requirement. 

Once  W  has  been  broken  down  into  its  many  component 
matrices  (4.1.4),  the  computational  procedure  becomes  much 
shorter;  in  fact  the  computing  time  required  is  g(log*q)/q2 
of  the  time  required  when  using  the  DFT  on  the  full  data 
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set . 

All  computations  are  done  using  complex  functions  such 

that,  for  real-valued  xK,  the  cosine  amplitude  is 

5-i 

ReC^j)  -  £  C  0  s  (  2  TV  K  j/s  )  (4.1.5) 

K *  o 

and  the  sine  amplitude  is 

9-1  N, 

=“£  Sl“  (ZTr*j/3)  .  («.1.6) 

It  is  easily  seen  that  for  frequency  zero  and  the 
Nyguist  frequency ,  the  imaginary  part  of  the  amplitude  will 
be  zero  when  using  real  data  and  an  FIT  technique  based  on  a 
power  of  two.  Also*  the  real  parr,  of  the  amplitude  is 
symmetric  about  the  Nyguist  frequency*  f  M  =  *  such  that. 

Ee  =  Re  ly^r)  for  r=i  ,2, . . . .  ,4-1  * 

and  the  imaginary  part  is  anti-symmetric  about  the  Nyguist 
frequency  such  that 

Im(y^_r)  =  ~XmiyK<r)  for  r  =  1 *  2,  .  . .  .  1 . 

The  values  for  the  cospectrum  are  calculated  using 

P;  =  2  { (Ee  t y j  ) )  2  +  (HMY;))2}  *  (4.1.7) 

Since  tne  data  are  standardized  by  removing  the  mean 
and  dividing  by  the  standard  deviation  prior  to  analysis, 
the  square  of  the  amplitude  is  recovered  from  (4.1.7)  by 
using 

=  2  Pj  " I  (4.1.8) 

0 


where  c-2  is  the  variance  of  the  g  data  values  and  <xr2-1  is 
the  unit  variance  of  the  standardized  values. 
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If,  as 
the  data  set 
factor  must 
added;  thus 


is  frequently  the  case,  the  number  of  values  in 
is  not  a  power  of  two,  a  second  multiplication 
be  used  in  order  to  compensate  for  the  zeroes 


(4.1.9) 


where  N  is  the  number  of  values  in  the  original  data  set. 
Therefore*  we  have  that 

A  J  2Pj0'13 

Aj  '  - H -  (4.1.10a) 

or  $  =  2P,a;*  (4.1.10  b) 

where  is  the  variance  of  the  original  N  data  values. 

The  corresponding  phase  and  frequency  are  readily 
obtained  by  using 


tan 


X  r*  (  d’t. S  ) 
R  e  { Afj) 


(4.1.11) 


and 


4 

J 

$  At 


(4.1.  12) 


where  At 

is  the 

time 

interval 

between 

data 

values  and 

3  0,1,.«.. 

,q/2.  The 

phase 

angle  <$j 

refers  to 

the 

phase  lag 

of  a  cosine  wave. 


4.2  Filters 

Two  filters  were  programmed  as  options  in  calculating 
the  FF T  of  tne  data.  The  first  was  a  tapering  filter 
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(Subroutine  TAPER)  which  would  taper  the  ends  of  the  data 
sample  in  order  to  avoid  pronounced  discontinuities  at  the 
ends-  The  second  was  a  Daniell  filter  (Subroutine  DANIj 
which  would  smooth  the  output  values  for  easier 
interpretation.  The  theory  for  both  filters  is  described  by 
Kana-sewich  (1^75). 

The  tapering  filter  consists  of  a  pair  of  cosine  bells. 
The  weights  for  one  half  of  the  filter  are  calculated  using 

Vi  =  T  (  I  +  co  s  jrJj^pL )_)  (4.2.1) 

with  iii-G,  1 1 . .  . .  fLf  and  are  applied  using 

l/m  and  ~  V-  a*-™ 

where  N  is  the  total  number  of  values  and  m=0/ 1 , • . . . #L . 

The  shape  of  the  filter  is  given  m  Figure  4.1. 


Figure  4.1.  Taper  filter. 

The  Daniell  filter  has  a  rectangular  or  box  car  shape 
response  in  tne  frequency  domain  as  shown  in  Figure  4.2. 

Kanasewich  recommends  using  the  Daniell  filter  on  the 
periodogram  calculated  by  the  FFT  since  it  reduces  the 
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Figure  4.2.  The  Daniell  or  rectangular  spectral  window. 

(From  Kanasewich,  1975. j 


variance  of  tne  periodogram. 

In  order  that  g  be  a  power  of  two,  the  data  must  be 
augmented  by  zeroes;  this  results  in  a  smoothing  of  the 
spectral  window  as  shown  in  Figure  4.3  for  N/g=0.7  and  1.0, 
with  r=7.  The  shape  of  the  window  is  given  by 


/ 

2  r+  i 


sin *  (  M to  - 
2  rrr  5<K)Z  (  LO  -  *  tj3j  J 


(4.2. 2) 


The  spectral  estimate  now  becomes 

?( U)  -  i  F(u-j)  F'(U-j)  (4.2.3) 

where  F  is  the  F FT  or  the  data,  F*"  is  the  complex  conjugate 
of  F  and  U=0, 1 g/2 .  The  factor  1/(2r+1)  normalizes  the 
window . 

Special  consideration  must  be  given  to  the  weighting 
factors  when  the  filter  is  applied  at  or  near  the  endpoints 
since  the  required  number  of  values  of  the  periodogram  is 


not  available  to  cover  the  full  width  of  the  filter.  In 
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Figure  4.3.  Danieli  window  for  N/g=0.7  and  1.0 ,  and  r=7. 

iFroin  Kanasewich,  1975.) 

order  to  avoid  this  problem,  the  end  points  were  ignored  and 
the  filter  was  applied  only  as  needed  to  obtain  the  required 
number  of  values  of  the  periodogram. 

4.3  Programs  and  Data  Manipulation 

The  prime  subroutine  used  in  tne  FFT  analysis  was 
FASTFT.  The  auxiliary  subroutines  required  were  TAPER, 
VARSD,  NLOGN,  POLAR,  DANI  and  DANY ;  subroutine  VARSD  will  be 
discussed  in  chapter  6.  All  programs  are  listed  in  Appendix 
B. 

Prior  to  the  FFT  analysis,  the  data  were  manipulated  as 


follows: 
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1-  The  data  mean  was  removed. 

2-  If  desired,  the  data  were  tapered  to  zero  at 
both  ends  using  a  cosine  bell  filter. 

3-  Zeroes  were  added  at  the  beginning  and/or  the 
end  of  the  data  set  until  the  desired  size  of  2" 
was  achieved. 

4-  The  mean,  variance  and  standard  deviation  were 
calculated  for  all  2*  values  using  VARSD. 

5-  The  data  were  standardized  using  the 
relationship 

Y\  =  (X;-EM)/SD  (4.3.1) 

where  EM  is  the  mean,  SD  is  the  standard 
deviation,  and  i=1,2,....,2n. 

6-  The  resultant  values  of  yL  were  assigned  to  the 
real  part  of  a  complex  array  X. 

The  actual  calculations  of  the  EFT  were  performed  by 
the  subroutine  NLOGN.  This  subroutine  appeared  in 
Robinson  (196  7)  « 

NLOGN  calculates  the  EFT  of  the  data  and  assigns  the 
results  to  the  complex  array  X.  This  subroutine  is  capable 
of  calculating  the  reverse  operation  rrom  transform  to  real 
data;  this  is  accomplished  by  specifying  the  index  of  the 
exponential  function  to  be  +1.0.  In  calculating  the 
transform,  the  index  is  specified  as  -1.0. 

The  complex  array  was  brolen  down  into  its  real  and 
imaginary  parts  from  which  the  power  spectrum  and  phase  were 
obtained  using  the  subroutine  POLAR.  The  subroutine  POLAR 
was  also  obtained  from  Robinson  (1967). 

The  square  of  the  amplitude  was  calculated  for  each 
value  in  the  power  spectrum  using  (4.1.10) . 

If  the  results  obtained  at  this  step  did  not  require 
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filtering,  they  were  printed  as  shown  in  Table  B. 15  in 
Appendix  B. 

The  rear  iCOS)  and  imaginary  (SIN)  parts  of  the  power 
spectrum  and  the  complete  power  spectrum  could  be  filtered 
using  a  Daniell  filter  ^Subroutine  DANI) ,  as  described  in 
the  previous  section.  The  resultant  smoothed  values  of  COS 
and  SIN  were  used  primarily  to  calculate  the  phase  values. 
The  frequency  was  filtered  in  a  slightly  different  manner  by 
using  a  modifred  Daniell  filter  ^Subroutine  DANYj  ;  the  power 
spectrum  was  used  as  a  weighting  factor  in  order  to  obtain  a 
better  estimate  of  the  true  peak  frequencies.  Both 
subroutines  nave  a  variable  filter  size  as  explained  by  the 
comment  cards  in  the  program  listing. 

All  values  obtained  were  printed  [see  Tables  B.16  and 


B.17  in  Appendix  B) . 


CHAPTER  5 


MAXIMUM  ENTROPY  METHOD 

5.  1  Theory 

J.P.  Burg  (1967,  1968)  proposed  a  new  spectral  analysis 
technique  called  the  Maximum  Entropy  Method  (MEM) .  The  first 
publication  discussing  MEM  appeared  in  1967,  at  the  same 
time  as  the  publication  promoting  the  EFT  was  released. 

This  technique  was  subsequently  extended  by  other 
authors,  in  particular  by  Parzen,  1969  and  Lacoss,  1971. 
T.J.  Ulrych  (1972)  gave  considerable  impetus  to  MEM  by 
obtaining  amazingly  good  results  with  che  technique.  Some  of 
these  results  were  later  challenged  by  Chen  and 
Stegen  (1974).  However  MEM  does  have  desirable  properties; 
the  positive  and  negative  aspects  of  this  method  will  be 
examined  in  Chapters  7  and  8. 

In  very  recent  time,  the  use  of  MeM  has  found  its  way 
into  meteorological  research  (Barn,  1976;  Dyer,  1976; 
Mason,  19  76)  . 

As  often  happens,  the  emergence  of  something  new  is 
soon  followed  by  variations  on  the  main  theme.  Ulrych  and 
Bishop  (1975)  provide  two  methods  lor  calculating  the 
spectral  estimates  using  MEM,  and  Andersen  (197U)  presents  a 
flow  chart  for  one  of  these  two  methods. 
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Kanasewich  (1975)  discusses  the  method  presented  by 
Andersen  (197  4)  and,  as  indicated  in  a  personal 
conversation,  recommends  this  approach. 

MEh  is  a  recursive  method  which  bears  greater 
resemblance  to  the  autocovariance  (ACV)  method  of  spectral 
analysis  than  to  the  FFT  method.  The  following  description 
of  MEK  is  taken  from  Chen  and  Slegen  (1974),  and  Andersen 


(1974)  . 

The  entropy  of  a  Gaussian  band-limited  time  series  is 
proportional  to 


(5. 1.  1) 


where  P(f)  is  the  power  spectrum  and  fw  is  the  Nyguist 
frequency.  Based  on  the  idea  that  the  spectral  estimate  must 
be  the  most  random  (i.e.  maximum  entropy)  of  any  power 
spectrum  that  is  consistent  with  the  measured  data. 
Burg  (1967,  1968)  maximized  the  entropy  shown  in  (5.1.1) 
through  the  use  of  Lagrange  multipliers  2ln  under  the 
constraints  of 

a  f  -N<n<N  (5.1.2) 

where  z  =  exp(i2rrf  at)  ,  At  is  the  sampling  interval,  and  <j>n 
is  the  autocorrelation  with  time  lag  naf.  Thus 


N 


a 


lf*0 


(5.  1.3) 


which  in  turn  yields 


' 


where  the  Lagrange  multipliers 
satisfying  the  constraints 


(5.1.4) 


are 

to 

be  determined 

by 

(5.  1 

.2)  . 

Now  „  (5.1.2) 

is 

equivalent  to 


(5.1.5) 


Since  P  (f )  is  real  and  nonnegative  and  f  tl  - ^at  * 
can  be  rewritten  as 


(5.1.4) 


_ Pn  _ 

gX  -  izrr  *i  &  t 

^  2—  ^  ^ 
n  ~l 


(5.  1.6) 


where  M  replaces  K  and  is  the  lag  value  chosen  to  optimize 
the  results. 


MEM  consists  of  calculating  a  set  of  filter 
coefficients  amn  and  a  constant  .  From  these  values  the 
power  spectrum  is  estimated  by  (5.1.6)  with  f  being  limited 
to  the  Nyquisr  interval f 
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so  that  Pm  is  the  output  power  of  the  m+1  long 
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.  Equation 

iteratively  by  stepwise  increase  of  the 
from  £)~1  to  m.  For  m=0,  Pe  is  estimatea  by 


(5. 1*7) 

matrix 


is  solved 
dimension 


(5®  1 «  8) 


In  general,  after  solution  of  (5.1.7)  for  m-1,  the  next 
step  m  involves  the  determination  of  a  set  of  m+2  unknowns 
(a*,,  ,  . .  «  .  ,  aWh);  <f^ ;  P*,)  from  the  m+  1  equations. 

According  to  Burg  (1968)  the  following  conditions 
should  be  used:  the  average  output  power  of  the  m* 1 
prediction  error  filter 


_  /  r~  *vj  rv)  jf""] 


is  minimized  with  respect  to  the  single  parameter  a*,*,.  The 
dependence  of  the  other  filter  coefficients  on  a*,™  is 
determined  by  the  m  lower  equations  in  (5.1.7)  having  the 
solutions 

-/  K  ~  ^  YA  r*  '  ^  to-l  n*'K  (5.1.10) 

for  k=1 , 2 m- 1  and  using  am0=~1  and  aMK=0  for  k>m. 

Using  (5.1.10),  (5.1.9)  can  be  written  as 
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2  (W-wi) 


N-M  f  \  x  /  i  \  2 

X.  L  (  ^  m  ^ yn  t)  ^  (  b Mi  _ 
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(5.1.11) 


where 


w>  v* 

~  am-lyA-K  't'i  +  m-K  /  (5.1.12a) 
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(5.  1. 12b) 
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Since  and  b1*,*.  are  independent  of  a^m/ 


—  TTm  =  o 


(5 .  1 . 13) 


Thus  Kn  is  a  minimum  with  respect  to  this  yields 
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Useful 

recursion 

formulas 

for  the 

arrays 

bmt  and  b£t 

derived 

bj;  means  of 

(5.1.1  0) 

and  (5.  1 

.12): 

hm  t  ~  brii  •£ 

^  ft*  -1  VH  -  1 

b  yyi  -i  t 

/ 

(5. 1. 15a) 

1  ^  M  -l  m  -1 

b  m-l  t+l 

* 

( 5 .  1  .  15b) 

It  is  seen  that  the  arrays  bmt.  and  b/L,*.  are  constructed 
the  arrays  bm_, e  and  b*n„ft  by  a  simple  linear  operation, 
starting  values  are 


hot  =  b't  =  x,, 

for  t=1tf 2, . . . . , N.  But  since  these  arrays  are  not  used  in  the 
calculation  procedure,  the  values  for  m= 1  are  used  instead 

b>t  =  xt  and  b*t  =  xt  +  ,  .  (5.1.16} 

By  inserting  (5.1.10)  into  (5.1.7),  the  recursive 
formula  for  P ^  is  found  to  be 

~  Pm  - 1  ( ^  “  3-wiy^)  »  (5.1.17) 


5.2  AJcaike  Filial  Prediction  Error 

The  number  M  of  filter  coefiicients  required  to 
optimize  the  results  obtained  from  MEM  can  be  determined  by 


two  methods. 
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One  metnod  is  based  on  a  window-closing  technique 
similar  to  the  one  described  by  Jenkins  and  Watts  (1968)  for 
the  variance  spectrum  analysis  method. 

A  superior  method  has  been  devised  by  Akaike  (1969a , 
1969b,  1970).  This  method,  now  called  the  Akaike  Final 
Prediction  Error  (FPE) ,  has  been  suggested  for  use  with  the 
AC V  method  (Akaike,  1971)  and  with  multivariate  time  series 
(Fryer  et  al.,  1975),  as  well  as  for  MEM. 

Since  the  emergence  of  the  Akaike  FPE,  many  authors 
have  incorporated  it  into  the  MEM  (Akaike,  1969a,  1969b, 
1970;  Gersh  and  Sharpe,  1973;  Jones,  1973;  Ulrych  and 
Bishop,  1975) . 

The  calculation  of  the  Akaike  FPE  is  readily  obtained 
by  using 


(5.2.  1} 


where  N  is  the  number  of  data  values  analysed,  M  is  the 
number  of  filter  coefficients,  and  S<~  is  the  residual  sum  of 
squares  given  by 


(5.2.2) 


in  which  aw^  are  the  filter  coefficients,  are  the  data 
values  with  tne  data  mean  removed,  and  xM=0  for  n<0. 

The  objective  is  to  fit  the  autoregressive  model  of 
order  M  by  the  least- sg uares  method;  'this  requires  the  mean 
square  of  residuals  (5.2.2)  to  be  minimized  with  respect  to 
{amM; m= 1 , 2, . . . . , Mj .  It  is  generally  observed  that  S& 
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decreases  and  (N+M+ 1) /  (N-M- 1 )  increases  as  the  value  of  M  is 
increased.  Thus  the  FPE  tends  to  be  too  large  when  a  large 
value  of  M  is  used.  The  minimum  value  of  the  FPE  provides  a 
balance  between  the  order  of  the  autoregressive  model  and 
the  size  of  tne  mean-square  prediction  error. 

An  estimate  (FPF)M  of  the  autoregressive  model  is 
calculated  by  (5.2.1) .  An  initial  value  is  calculated  using 


(5.2.3) 


where  P0  is  given  by  (5.1.8). 

For  the  purpose  of  comparison  of  the  magnitude  of 
(FPE)m  with  (1  PE)e  ,  the  relative  value 


(5.2.  4) 


can  be  used. 

5-5  Programs  and  Data  Manipulation 

Two  subroutines  were  used  to  calculate  the  power 
spectrum  using  MEM;  they  are  BURGME  and  SMSQR.  Both  are 
listed  in  Appendix  B. 

The  data  were  standardized  before  being  passed  on  to 
the  subroutine  BURGME  using  ( 4 .3.1)  for  i  1,2,.«.«,N,  w  h  e  i_  e 
N  is  the  totar  number  of  data  values. 

The  subroutine  BURGME  combines  the  theory  of  the  MEM 
with  the  Akaike  FPE.  The  theory  and  the  flow  chart, 
presented  in  Andersen's  paper,  were  used  extensively  in 
writing  the  subroutine. 


The  Akaike  FPE  has  been  incorporated  as  an  option  such 
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that  the  minimum  value  of  (FPE)^  can  ne  automatically  found 
and  the  associated  filter  coefficients  retained.  However, 
the  number  of  filter  coefficients  can  be  specified, 
bypassing  the  calculation  of  the  FPE. 

An  expansion  factor  has  also  been  incorporated  into  the 
program  in  order  to  expand  the  lower  end  of  the  spectrum  for 
better  definition  of  peats.  If  Q  is  this  expansion  factor, 
then  the  printout  covers  the  range  of  frequency  from  f =0  to 
i=f w/Q  where  i #  is  the  Nyquist  frequency.  In  this  manner, 
there  is  no  aliasing  of  values  in  the  range  f =f */QF . . . „ , fn 
onto  the  range  f=0,.„..,f u/Qi  aliasing  of  frequencies  beyond 
fw  must  nevertheless  be  investigated;  this  will  be  done  in 
Chapters  7  and  8. 


The  subroutine 

SdSQR 

calculates  the 

residual  sum 

of 

squares  S*  given  by 

( 5  *  2  e  2 ) 

and  required 

by  (5.2.1) 

if 

Akaike's  FPE  is  used  to  obtain  the  optimum  number  of  filter 
coefficients . 

Sample  print-outs  are  shown  in  Tables  B.18  and  B. 19  in 


Appendix  B. 


CHAPTER  6 


AUXILIARY  PROGRAMS 

6 , 1  Main  and  Plotting 

The  Earn  program  (MAIN)  linns  together  all  the 
subroutines  mentioned  in  Chapters  4  and  5.  It  also  allows 
for  plotting  the  raw  data  and  the  power  spectra  from  both 
the  PFI  and  MEM  analyses. 

The  plotting  is  accomplished  with  subroutine  PLTG  which 
is  a  slightly  modified  version  of  a  program  available  in  the 
Computing  Services  Program  Library  of  the  University  of 
Alberta,  The  subroutine  PLTG  is  not  included  in  Appendix  B 
as  it  is  too  lengthy  and  could  prove  difficult  to  implement 
on  other  computer  systems. 

The  subroutines  PLTDAT,  PLTFFT  and  PL'TMEM  provide  an 
interface  between  MAIN,  FASTFT,  BURGEE,  and  PLTG.  These 
three  interface  subroutines  isolate  most  of  the  commands 
required  when  producing  a  plot  and  provide  an  easy  method 
for  inserting  a  different  plotting  subroutine. 

If  no  plotting  facilities  are  available,  the  program 
need  only  be  slightly  altered  by  deleting  the  subroutines 
PLTDAT,  PLTFFT  and  PLTMSM  and  by  removing  all  CALL 
statements  to  these  three  subroutines  rrom  MAIN.  The  three 
subroutines  are  listed  in  Appendix  B. 
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6 . 2  Mean  and  Variance 

The  subroutine  VABSD  is  used  to  calculate  the  mean  ^ a }  * 
variance  (cr2 j  and  standard  deviation  (<r)  of  the  data.  It  is 
called  once  in  MAIN  in  order  to  obtain  these  three  values 
from  the  raai  data.  It  is  also  called  in  FASTFT  in  order  to 
obtain  a  new  set  of  values  from  the  expanded  data  set. 

The  mean  is  calculated  using 


and  the  variance  is  calculated  using 


(6. 2. 2) 


Statistically s  the  values^  and  cr2  are  known  as  biased 
estimators  of  the  true  population  values. 

A  listing  of  VAfiSD  can  be  seen  in  Appendix  3. 

6 • 3  Notch  filter 

It  was  anticipated  that  some  method  would  be  required 
to  eliminate  the  very  strong  annual  cycle  when  analysing  a 
full  set  of  temperature  data  at  one-month  intervals. 

The  filter  designed  for  this  purpose  is  called  a  notch 
filter;  if  is  also  referred  to  as  a  rejection  filter 
(Kanasewich*  1975}  or  a  recursion  filter  in  the  z-plane 
(Shanks*  1967} . 

The  following  theory*  as  described  by  Kanasewich  (1975) 
and  Shanks  (1967)  *  shows  how  the  notch  filter  is  used  to 
remove  the  annual  cycle.  The  program*  as  written  (see 


Appendix  B) ,  can  be  used  to  eliminate  any  frequency. 

Consider  a  unit  circle  in  the  imaginary  plane.  The 
positive  imaginary  plane  progresses  from  f=0  to  £=f the 
Nyquist  frequency. 

The  general  form  of  the  filter  equation  is  given  by 

Fh?)  -  -L - — - -  (6.3.1) 

(  2  -  S,)  (  *  -  ?*) 

where  z,  and  z*  are  the  zeroes  and  z3  and  zH  are  the  poles 
of  the  filter. 

The  zeroes ,  or  discontinuities,  of  the  filter  occur  at 

z  =  cos  9  ±  i  sin  © 

where  Q  is  the  position  on  the  unit  circle  of  the  unwanted 
frequency.  Thus#  for  £«,  the  frequency  of  the  annual  cycle, 
the  value  of  0  is  determined  by 

©=  JL  .  IS0°  -  -th-  •  I  80°  =  30° 

h  • 5 

assuming  the  data  interval  to  be  one  month. 

The  zeroes  of  the  filter  occur  at 

z,  —  0.8660  +  10.5000  (6.3.2a) 

and  z2  =  0.  8660  -  10.5000  .  (6.3.2b) 

However,  when  the  poles  of  the  filter  coincide  with  the 
zeroes  on  the  unit  circle,  the  response  is  very  poor.  It  is 
necessary  to  place  the  poles  near,  but  not  on,  the  unit 
circle,  say  az  | z {  =  1.001.  Thus  the  poles  are  located  at 

z3  =  0.8669  +  i0.o005  (6.3.3a) 

and  Zv  =  0.8669  -  i0.5005  .  (6.3.3b) 


Substit  uuing  (6.3.2)  and  (6.3. a)  info  (6.3.1),  the 
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filter  equation  becomes 


F(z)  =  0.'>980(i'-l.73ZOi1])m  (6.3.4, 

/  -  /.  7  3  O  3  z  +0.99802* 

which  incorporates  a  static  gain  factor  to  insure  a  gain  of 
unity  at  the  Nyquist  frequency. 

Next,  the  filter  is  applied  to  the  data  using 

Y  (z)  =  F  (z)  X(z)  (6.3.5) 

where  X  are  the  raw  data  and  Y  are  the  filtered  data.  The 
recursion  equation  is  obtained  by  substituting  (6.3.4)  into 
(6. 3c  5),  by  multiplying  both  sides  of  (6.3.5)  by  the 
denominator  and  by  rearranging  to  obtarn 
Y  (zj  =  0.99  80  (1-u  7320z  +  z2)X  (z } 

+  zY  (z)  (1. 7303-0. 9980z)  .  (6.3.6) 

Whenever  z  occurs  in  (6.3.6),  the  input  or  output  is 
delayed  by  one  unit  and  similarly  z2  causes  a  delay  of  two 
units.  Thus,  the  recursion  equation  is 

y„  =  0.9980  (  Xy,  -  1.732  0XV,..,  +  x*-*  ) 

+  1 . 7  3G3yv1_,  -  0.9980y„.z  (6.3.7) 

where  xM  is  the  data  input  and  yn  is  the  output  from  the 
filter. 

The  amplitude  response  of  this  filter  is  shown  in 
Figure  6.1  in  which  120  frequency  values  are  assumed  to  lie 
in  the  interval  0<f<f w * 

The  recursion  equation  (6.3.7)  is  programmed  to  filter 
the  data  in  a  cascade  fashion.  A  cascaae  consists  of  passing 
the  filter  twice  over  the  data  set,  once  in  a  forward 
direction  from  x, 


to  xw  and  once  in  a  backward  direction 


/ 
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Figure  6.1c  Amplitude  response  of  the  notch  .filter 

for  the  annual  cycle. 


from  xw  to  x,  *  The  cascade 
data.  Also,,  a  sufficient 
the  ends  of  the  input  data 
initial  pass  of  the  filter 
The  program  .  for  the 
spectral  program. 


prevents  a  phase  shift  in  the 
number  of  zeroes  must  be  added  to 
to  allow  the  output  from  the 
to  become  significantly  small, 
notch  filter  is  separate  from  the 
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CHAPTER  7 


RESULTS 


7 . 1  Program  P erf ormance 

Part  of  the  development  of  the  programs  involved  tests 
which  were  performed  using  synthetic  data.  These  data 
consisted  of  different  combination: 
periodicities.  Every  generated 


tions  of 

f  ive 

to  seven  exact 

ated  wa 

ve 

had  a  specific 

angle,  i 

'h  e 

waves  had  small 

es  and 

larg 

er  amplitudes  at 

iods  ranged 

from  i 00  to  2 

Hyears,!«  The  length  of  data  used  was  equivalent  to 
approximately  70  years  of  real  data.  Thus,  when  using  the 
modified  periodogram ,  zeroes  were  added  to  the  generated 
data  as  if  real  data  of  finite  length  were  being  used.  These 
generated  data  were  analysed  using  the  modified  periodogram 
method  and  MEd. 

The  modified  periodogram  performed  quite  well  and  was 
able  to  analyze  the  periodicities  provided.  Periods  greater 
than  approximately  one-third  the  length  of  data  were  not 
resolved  as  accurately  as  the  shorter  periods.  Periods 
approximately  equal  to  or  greater  rhan  the  length  of  data 
were  detected  at  the  first  frequency  value  above  zero,  but 
no  values  of  frequency,  amplitude  or  phase  could  be 
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determined  with  precision.  The  values  of  the  spectral  peaks 
exceeded  the  99.9%  significance  level  for  the  resolvable 
waves  present  in  the  synthetic  data  (srgnif icance  levels  are 
discussed  in  the  next  section).  The  spectral  peaks  were 
analyzed  at  tne  output  frequencies  which  were  closest  to  the 
frequencies  specified  in  generating  the  data.  The  weighted 
freguency  proved  useful  in  this  regard  because  it  provided 
interpolated  values  of  frequencies  which  were  better 
estimates  of  the  actual  values.  Similarly {  the  smoothed 
phase,  obtained  using  a  three-point  Daniell  filter  on  the 
real  and  imaginary  parts  of  the  Fourier  transforms,  also 
provided  better  estimates  of  the  phase  angles  than  the 
unsmoothed  phase  angles. 

Tapering  was  performed  on  the  data  in  order  to  assess 
the  effects  of  truncating  the  data  set.  No  improvements  were 
observed  in  the  resolution  properties  of  the  per iodogram. 
All  periods  within  the  range  of  detectable  frequencies  were 
revealed  wit a  equal  or  slightly  superior  results  when  no 
tapering  was  applied  to  the  ends  of  the  data. 

The  Maximum  Entropy  Method  did  not  prove  to  be  as 
accurate  or  as  convenient  as  the  modiried  per  iodogram.,  Many 
observations  made  about  MEM  in  this  study  were  also 
discussed  by  Chen  and  Stegen  (1974). 

The  spectral  output  of  MEM  was  more  difficult  to 
interpret.  It  required  considerable  time  and  care  when 
determining  che  power  resolved  for  a  certain  freguency. 


since  it  was  not  the  value  of  the  pear  which  provided  an 
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estimate  of  the  power  bat  rather  the  area  under  the  curve, 
for  example*  the  maximum  single  value  of  variance  attributed 
to  two  periods  of  equal  amplitude*  among  a  combination  of 
five  generated  periods*  differed  by  a  factor  of  twenty  in 
the  same  analysis.  Thus,  confidence  limits  could  not  be 
applied  as  easily  to  MEM  results  as  was  possible  with  the 
modified  periodograra  results.  At  this  time*  the  true  meaning 
of  the  numerical  values  calculated  by  MEM  is  unknown.  The 
value  of  a  spectral  peak  is  subject  to  wild  fluctuations 
when  slightly  different  data  sets  are  used.  The  periodogram 
method  was  not  subject  to  such  wild  fluctuations  and  the 
maximum  single  value  of  variance  attributed  to  periods  of 
-equal  amplitude  differed  by  less  than  a  factor  of  two.  MEM* 
unlike  the  periodograra  method,  must  be  calibrated  using  test 
data  before  being  used  on  real  data. 

The  resolution  properties  of  MEM  were  highly  dependent- 
on  the  number  of  filter  coefficients  used,  especially  if  MEM 
was  used  to  determine  the  value  of  a  very  low  frequency.  The 
number  of  coefficients  required  was  usually  higher  than  the 
number  determined  by  the  Akaike  FPE.  Thus,  MEM  required  a 
process  similar  to  the  determination  of  the  optimum  lag 
value  for  the  ACV  method. 

Also,  the  spectral  peaks  analyzed  were  not  necessarily 
at  the  proper  frequency.  Chen  and  Stegen  (13?a)  found 
frequency  shirts  which  resulted  from  a  combination  of  the 
initial  phase  and  the  data  length.  At  this  stage,  MEM  can 
only  provide  an  indication  of  low  frequencies  present  in  the 
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data.  However,  a  recent  paper  by  Ulrych  and  Clayton  (1976) 
discusses  a  method  by  which  the  frequency  shift  may  be 
reduced  considerably .  This  paper  was  brought  to  the 
attention  of  the  author  by  Dr.  T.J.  Uirych  during  a  private 
conversation.  Unfortunately,,  this  conversation  took  place  at 
a  very  late  stage  of  this  study  and  no  time  was  available  to 
make  the  modifications  necessary  to  improve  the  method. 

MEW  was  round  to  be  quite  costly  in  computer  time,, 
especially  when  using  long  data  sets.  The  time  required  when 
calculating  the  Akaike  FPE  was  also  significant  and  was 
mainly  due  to  the  calculation  of  5^  using  (5.2.2).  A 
modification,  obtained  more  by  accident  than  by  design,  was 
discussed  witn  Dr.  Ulrych.  He  outlined  the  justification  for 
replacing  the  value  of  obtained  from  (5.2.2)  by  the 
latest  value  of  Pp,  obtained  by  (5.1.17).  This  substitution 
is  justified  since  both  and  PM  are  measures  of  the 
variance  of  the  error  between  the  calculated  set  of  filter 
coefficients  and  the  ideal  set  of  filter  coefficients. 
Substituting  PM  for  the  calculation  of  would  reduce  the 
computational  time  required. 

The  notch  filter  was  found  to  be  quite  effective  in 
reducing  or  eliminating  an  unwanted  frequency.  A  "before  and 
after"  example  of  the  effect  of  the  filter  is  shown  in 
figures  7.1  and  7.2.  The  upward  shift  of  the  plot  in  Figure 
7.2  results  from  the  redistribution  of  that  portion  of  the 
normalized  power  which  has  been  released  by  eliminating  the 
annual  cycle. 
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Figure  7.1.  Periodogram  analysis  of  FV R  monthly  temperatures 

with  no  filtering  appli.ed  to  the  annual  cycle. 


Figure  7.2.  Periodogram  analysis  of  FV*  monthly  temperatures 

with  the  annual  cycle  removed. 
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7 • 2  Data  Analysis  and  Evaluation  Procedure 

The  results  of  the  climatological  analysis  are 
discussed  under  four  main  headings  *  This  section  deals  with 
the  procedure  used  to  analyze  and  evaluate  the  results. 
Section  7.3  deals  with  results  affecting  Alberta  as  a  whole 
based  on  air  stations  analyzed.  The  remaining  sections  of 
this  chapter  discuss  each  station  on  a  regional  basis. 

The  primary  method  used  was  the  periodogram  analysis. 
Unf ortunately e  MEM  was  used  very  little  for  reasons  which 
have  been  detailed  in  the  previous  section. 

The  procedure  used  to  analyze  the  data  consisted  of 
performing  a  periodogram  analysis  on  various  data  groupings 
as  outlined  in  Section  1.3.  The  power  spectra  obtained  were 
plotted  on  a  Calcomp  plotter.  Both  smoothed  and  unsmoothed 
spectra  were  calculated.  The  smoothing  was  normally 
performed  with  a  three-point  Danieil  filter.  Nearly  600 
graphs  of  power  spectra  were  eventually  produced. 

Confidence  intervals  were  calculated  using  the 
assumption  of  a  white-noise  spectrum  and  a  Chi-square 
distribution.  The  value  of  the  white-noise  variance  was 
dependent  on  the  number  of  spectral  lines  analyzed  and  was 
easy  to  determine.  The  confidence  limits  were  obtained  by 
multiplying  the  white-noise  value  by  X ,_<&//  and  X^/y  tor 
the  upper  and  lower  limits,  respectively ,  where  y  is  the 
appropriate  number  of  degrees  of  freedom  and  oc  is  the 
probability  level  selected.  The  degrees  of  freedom  for  the 
periodogram  analysis  were  determined  by  using 
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where  N  is  the  number  of  independent  data  values  used  in  the 
analysis,  2M_1  is  the  number  of  spectral  lines  obtained  when 
the  N  data  values  are  supplemented  by  2M-N  zeroes  (see 
Section  4.3),  and  D  is  the  width  of  the  Daniell  filter  used 
on  the  periodogram  (D=1  when  the  periodogram  is  unsmoothed). 
The  number  of  data  values  is  reduced  by  two  degrees  of 
freedom  as  a  result  of  using  the  data  to  calculate  the  mean 
and  variance  values. 

The  degrees  of  freedom  for  the  unsmoothed  periodogram 
varied  between  one  and  two.  If  the  calculated  degrees  of 
freedom  were  1.5,  an  erroneous  result  could  have  been 
obtained  by  using  y  =  1  (over-restricting)  or  V-2  (under- 
restricting)  since  the  distribution  of  the  Chi-square 
function  varies  considerably  and  non-lmearly  between  these 
two  values  of  V.  For  this  reason,  grapns  were  plotted  of  the 
Chi-square  function  in  order  to  obtain  more  accurate 
estimates  of  the  confidence  intervals  for  both  unsmoothed 
and  smoothed  spectra. 

The  95%  and  5%  confidence  limits  and  the  white-noise 
spectrum  were  drawn  on  the  graphs  as  snown  in  Figure  7.3. 
All  spectral  lines  exceeding  these  limits  were  noted  along 
with  the  significance  level  reached.  Tne  levels  considered 
were  0.5%,  1*,  2.5%,  5%,  95%,  97.5%,  99%,  99.5%,  and  99.9%. 
The  spectrum  illustrated  in  Figure  7.3  is  typical  of  the 
spectra  obtained  in  this  study  and  it  demonstrates  the 
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Figure  7.3.  Periodogram  analysis  of  YXH  winter 
precipitation  showing  a  white-noise  spectrum  and 
a  90%  confidence  interval. 


validity  of  assuming  a  white-noise  model. 

The  analysis  for  all  months  and  for  each  month  of  the 
calendar  year  was  accomplished  using  rhe  data  as  collected. 
The  temperature  and  precipitation  data  were  averaged  and 
summed  respectively  o-ver  the  appropriate  time  interval  to 
obtain  the  annual  and  seasonal  values.  The  groupings  were 
made  as  follows: 

1-  Annual  -  January  to  December. 

2-  Winter  -  December,  January,  February. 

3-  Spring  -  March,  April,  May. 

4-  Summer  -  June,  July,  August. 

5-  Fall  -  September,  October,  November. 
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These  groupings,  as  opposed  to  the  full  sets  of  monthly 
values,  were  chosen  because  it.  was  felt  that  they  would 
eliminate  the  necessity  of  removing  the  annual  cycle  which 
generally  accounted  for  80%  or  more  of  the  data  variance. 
The  notch  filter  proved  capable  of  eliminating  this  period, 
however,  possible  side-effects  of  removing  such  a  tremendous 
portion  of  tne  variance  might  have  inrluenced  the  remaining 
results. 

Beaverlodge  was  chosen  for  an  in-depth  analysis  of  some 
of  the  results  obtained.  This  analysis  is  discussed  in 
Section  7.3. 

Periods  of  32  years  or  greater  were  considered  to  be 
very  doubtful  because,  in  most  cases,  only  64  spectral  lines 
were  analyzed.  An  example  of  such  an  analysis  is  shown  in 
Table  B.  15.  The  first  five  values  of  frequency  output 
correspond  to  periods  of  128,  64,  42.7,  and  32  years.  A 
significant  peak  with  a  period  of  32  years  cr  more  was 
viewed  as  possible  but  not  necessarily  definitive  evidence 
of  long  periodicities  present  in  the  data.  Such  significant 
periods  are  included  in  Section  7.3.  However,  a  discussion 
of  their  validity  is  reserved  for  Section  7.4. 

In  addition  it  is  often  difficult  to  decide  whether  a 
wave  is  signiiicant  by  itself  or  is  a  harmonic  of  some  other 
wave;  for  example,  significant  periods  of  8  and  4  years 
could  easilv  result  from  the  analysis  of  a  single  8— year 
wave  pattern  where  both  waves  are  needed  to  describe 
accurately  the  wave.  The  identification  of  such  a  harmonic 
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wave  would  simplify  the  task  of  finding  the  causes  of 
significant  periods. 

It  must  be  remembered  that  the  analyses  of  annual, 
seasonal,  and  individual  month  variations  for  BEA  and  MAN 
provided  only  32  spectral  estimates  between  zero  and  the 
Nygurst  frequency  while  64  estimates  were  obtained  for  all 
other  stations.  Also,  5%  of  these  estimates  were  expected  to 
exceed  the  95^  confidence  limit.  Therefore,  an  analysis  with 
three  estimates  exceeding  the  95%  significance  level  out  of 
64  spectral  estimates  did  not  necessarily  provide  meaningful 
results,  as  will  be  shown  in  the  remaining  sections  of  this 
chapter. 

7.3  Alberta  as  a  Unit 

A  province-wide  comparison  was  limited  mostly  to 
results  obtained  from  annual  and  seasonal  values  of 
temperature  and  precipitation.  This  limitation  was  due  to 
the  volume  of  results  generated  and  tne  time  restrictions 
involved.  A  more  detailed  analysis  or  individual  months  is 
discussed  in  the  last  three  sections  of  this  chapter. 

The  "spectral  gaps"  or  estimates  which  are  below  the  5% 
confidence  limit  are  not  discussed  at  this  time.  Instead, 
attention  is  focused  on  frequencies  or  groups  of  frequencies 
that  appear  to  account  for  relatively  large  amounts  of  total 
variance. 

Tables  7.1  to  7.5  show  the  results  obtained  for  the 
analysis  of  annual  and  seasonal  values.  Only  significant 
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Table  7.1.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  annual  values. 


Temper at  ure 

Precipit  at  ion 

p  i% 

)  99  97.5 

95 

99  97.5 

95 

FVR 

128.  64.0 

16.0 

14.2 

128. 

5.33 

BEA 

16.0^  64.0 

7.  11 

RAN 

none 

none 

yxd 

none 

10.7 

2.25 

LAC 

2.67 

2.61 

9.85 

9.14 

LTH 

2.67 

12.8 

4.27 

YXH 

64 .  O'4*  128. 

4.92 

18.3 

25.6 

14.2 

2.37 

MAN 

21.3 

64.0 

-*• 

indicates  that  the 

Deak 

r 

exceeded  the  99.5%  level 

• 

peaks 

(>95%)  ,  obtained 

from 

the  unsmoothed  periodogram 

analyses,  are  listed  in  these  tables.  The  reasons  for  not 
using  the  smoothed  periodogram  results  are  discussed  in 
Section  7.5. 

Few  results  for  any  one  particular  station  exceeded  the 
expected  number  of  significant  peaks.  One  such  case  is  in 
the  winter  precipitation  of  YXH  from  which  tour  significant 
peaks  were  resolved.  However,  periods  21.3  and  18. j  aj-e 
adjacent  peaks  (see  Figure  7.3)  which  are  likely  to  indicate 
a  single  period  of  approximately  19.7  years.  This  then 
reduces  the  number  of  peaks  to  the  expected  number  although 


Table  7c 2.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  winter  values. 


Temperat  ure 

Precipitation 

p  (%) 

99 

97.5 

95 

99  97.5  95 

FVR 

16.0 

128.  42.7 

64.0**  4.27 

BEA 

16.0 

3.20 

2.21 

none 

RAN 

3.28 

2.25 

2.37 

YXD 

3.28 

10.  7 

2.21 

3.20 

LAC 

none 

4. 74 

LTH 

6.  74 

3.28 

64.0 

YXH 

3.28 

2.46 

21.3+  2.78  128. 

18.3* 

WAN 

3.20 

16.0 

7.11 

3.39 

64 . 0* 

*  indicates 
**  indicates 

that  the 
that  the 

peak 

peak 

exceeded  the  99.5%  level, 
exceeded  the  99.9%  level. 

the  fact  remains  that  these  adjacent  peaks  are  both 
significant  at  the  99.5%  level  and  that  together  they 
account  for  26.67a  of  the  variance  of  the  data.  This  result 
is  considered  further  in  Section  7.6. 

Three  results  are  the  most  prominent  when  all  stations 
are  examined.  All  three  results  come  from  the  temperature 
analyses.  They  consist  of  a  3.2-  to  3.3-year  period  in  the 
winter  at  sn  stations,  a  3.8-  to  4-year  period  in  the 
spring  at  five  stations,  and  a  4.9-  to  5.1-year  period  in 
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Table  7.3.  Periods 

(years) 

of  significant  peaks  and 

significance 

levels 

(P,%)  for  spring  values. 

Temperat ure 

Precipitation 

P »%)  99  9  7.5 

95 

99  97.5  95 

FVR 

128. 

16.0 

14.2 

6.7a 

BEA 

16.0* 

non  e 

HAN 

3.88 

2.42 

128.  64.0 

YXD 

5.12 

4.00 

2.61 

none 

LAC 

3.88 

none 

LTH 

3.88 

none 

YXH 

123. 

25.6 

64. 0 

18.3 

4.00 

2.61 

6.74 

MAN 

2.46 

10.7 

*  indicates  that  the 

peak  exceeded  the  99.5%  level. 

the  fall  at  all  stations. 

The  ffiosi  prominent  result  is  the  approximate  five-year 
period  in  the  fall  temperature  values.  All  stations  have  a 
significant  peak  at  4.92  years  ana  five  of  these  stations 
have  a  second  significant  peak  at  an  adjacent  period  such 
that  the  true  period  can  be  thought  to  lie  somewhere  between 
the  two  peaks.  In  the  fall  values,  the  amplitudes  of  the 
significant  peaks  at  and  adjacent  to  the  4.92-year  peak 
account  for  10  to  30 %  of  the  variance  of  the  data  used.  This 
result  is  as  baffling  as  it  is  surprising  because  it 
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Table  7.4.  Periods  (years)  of  significant  peaks  and 
significance  levels  {P,%)  for  summer  values. 


P(%) 

Temperature 

Precipitation 

99 

97.5 

95 

99  97.5  95 

FVE 

32.  0 

128. 

8.00 

3.39 

42.7 

2.67 

BEA 

2.37 

16.0 

7.11 

2.13 

FAN 

42.7 

21.3 

non  e 

YXD 

128. ** 

32.  0 

10.7  12.8 

LAC 

25.  6 

3.80 

LTH 

12.8 

42.  7 

12.8 

YXH 

64 «  O'4"* 

128. 

2.3  7*  12.8 

42.  7 

12.  8 

4.00 

MAN 

none 

none 

+  indicates 

that  the  peak 

exceeded  the  99.5%  level. 

indicates 

that  the  peak 

exceeded  the  99.9%  level. 

provides  possible  support  for  the  results  obtained  by 
Georgiades  (1977)  while  having  no  known  physical  cause. 

This  five-year  period  was  investigated  as  to  possible 
sources  other  than  a  true  periodicity  present  in  the  data. 
This  investigation  was  carried  out  using  fall  temperatures 
for  BEA .  A  plot  of  the  fail  temperatures  is  shown  in  Figure 
7.4.  The  pericdogram  analysis  of  the  raw  data  is  shown  in 
Figure  7.5.  An  11-point  Tukey  filter  was  passed  over  the  raw 
data;  this  reduced  the  number  of  data  values  from  61  to  53. 
The  filtered  data  are  shown  in  Figure  7.6  and  the 
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Table  7.5.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P/%)  for  fall  values. 


Temperat  are 

Precipit  at  ion 

99  97.5  95 

99  97.5  95 

F  VS 

5.  12 

4.92 

non  e 

BEA 

4.92 

6.40 

2.56 

32.0  3.05 

SAN 

5.  12 

4.92 

2.84 

2.78 

YXD 

4 . 92* 

4.57 

LAC 

5.  12 

4. 92 

10.7  6.ao 

2.51 

LTH 

5.  12 

4.92 

7.53 

non  e 

YXII 

4.92* 

7.53 

16.0 

HAN 

4.92 

4.57 

9.14 

t 

indicates 

that  the 

peak  exceeded  the 

99.5%  level. 

periodogram  analysis  of  these  data  is  shown  in  Figure  7.7. 

Comparing  Figures  7.5  and  7.7  it  will  be  seen  that  the 
power  of  the  low  frequencies  was  reduced  and  the  power  of 
most  higher  frequencies  was  increased  slightly  as  a  result 
of  the  redistribution  of  the  total  power.  Thus  leakage  of 
power  from  the  low  frequencies  did  not  appear  to  contribute 
to  the  power  of  the  five-year  period. 

The  next,  problem  was  to  determine  if  aliasing  had  taken 
place.  The  power  spectra  of  daily  and  monthly  temperature 
values  were  not  included  in  this  thesis  due  to  the  size  o.t 
the  figure  required  for  adequate  display.  A  series  of  daily 
mean  temperature  values  was  readily  available  lor  BEA.  An 
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Figure  7,4.  BEA  fall  temperatures. 


Figure  7.5.  Periodograra  analysis  of  BEA  fall  temperatures. 
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YERR 

Figure  7.6.  Filtered  values  of  BEA  rail  temperatures. 


Figure  7.7.  Periodogram  analysis  oi  rhe  filtered  values 

of  BEA  fall  temperatures. 
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analysis  of  four  consecutive  years,  1950  to  1953,  showed 
that  daily  values  follow  a  red-noise  spectrum.  The  annual 
cycle  was  clearly  evident  and  no  other  frequency  was  found 
to  be  significant. 

The  analysis  of  monthly  values,  in  which  all  monthly 
mean  values  of  temperatures  were  used,  was  performed  after 
the  removal  of  the  annual  and  semi-annual  cycles  using  the 
notch  filter.  This  analysis  showed  significant  peaks 
associated  with  16. 7-,  1.82-,  0.30-,  and  0.208-year  periods. 
The  last  three  periods  are  beyond  the  Nyquist  frequency  of 
0.5  as  shown  in  Figure  7.5.  Using  a  linear  folding  method, 
the  three  periods  would  be  aliased  onto  the  2.22-,  3.28-, 
and  5.0-year  periods.  However,  the  variance  observed  at  each 
of  the  three  highest  frequencies  was  less  than  1%  of  the 
total  variance  of  the  data.  All  these  factors  are  discussed 
in  the  next  cnapter  before  drawing  conclusions. 

At  this  time,  a  few  more  facts  concerning  fall 
temperature  analyses  are  presented.  Of  the  eight  stations, 
three  records  (FVR,  LAC  and  LTH)  began  in  1908.  The  phase 
angles  computed  were  43°,  2°,  and  1°,  respectively.  The 
phase  angle  for  YXD,  beginning  in  1881,  was  51°,  and  for 
YXH ,  beginning  in  1884,  was  -82°;  using  an  approximate  value 
of  72°  for  the  phase  change  each  year  based  on  an  exact 
period  of  five  years,  the  phase  angle  for  YXH  in  1881  would 
have  been  62°.  The  phase  angles  obtained  by  adjusting  the 
values  for  3EA  and  RAN  from  1915  and  1905  to  1908  were  60° 
and  38°,  respectively.  In  order  to  compare  the  phase  angle 
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of  two  stations,  a  correction  of  72°  per  year  was  applied  to 
the  derived  pnase  shift.  This  correction  is  not  likely  to  be 
exactly  72°  such  that  30-  to  40-year  corrections  could  not 
easily  be  applied  in  order  to  compare  all  stations.  The 
similarity  of  the  resolved  phase  angles  is  nonetheless 
striking . 

The  amplitude  of  the  five-year  period  was  approximately 
1C  for  most  stations.  The  values  which  differed  the  most 
from  1C  were  0. 78C  at  YXD,  0.84C  at  YXH,  and  1.16C  at  MAN. 
The  data  variance  was  smallest  at  MAN  with  a  value  of  2.20C2 
such  that  this  period  accounted  for  30.4%  of  the  variance. 

The  values  obtained  from  the  analyses  of  individual 
months  also  provided  some  interesting  results  with  respect 
to  the  five-y*ear  period.  Significant  peaks  associated  with 
this  period  were  found  for  six  stations  in  the  September 
data,  and  for  five  stations  in  the  November  data,  but  no 
station  revealed  such  a  peak  in  the  October  data.  It  was 
found  upon  closer  examination  that,  for  all  stations,  the 
percentage  of  variance  explained  in  the  October  data  by  an 
approximate  five-year  period  was  equivalent  to,  or  less 
than,  the  percentage  explained  by  a  white-noise  spectrum. 

The  analyses  of  summer  values  revealed  no  single 
frequency  which  was  as  frequently  significant  as  the 
frequencies  revealed  for  the  other  seasonal  values. 
Generally,  tne  significant  peaks  occurred  at  very  low 
frequencies.  Such  a  result  is  very  difficult  to  interpret 
using  the  periodogram,  and,  as  discussed  in  Section  7.1,  the 
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MEM  analysis  would  add  little  information. 

Both  stations  which  did  not  exhibit  a  3.2-  to  3.3-year 
period  in  winter  temperatures  at  the  95%  confidence  level 
did  have  such  a  period  at  the  90%  confidence  level. 
Similarly  all  three  stations  which  did  not  exhibit  a  3.8-  to 
4-year  period  in  spring  temperatures  did  have  such  a  period 
at  the  85  to  90%  confidence  level.  Considerable  doubt  is 
cast  on  the  statistical  significance  of  these  results  by 
going  below  tne  95%  level.  However,  because  so  many  stations 
had  these  results  in  common,  the  presence  or  absence  of 
these  periods  at  the  other  stations  could  strengthen  or 
weaken  the  hypothesis  that  an  approximate  3.2-year  period  in 
winter  temperatures  and  an  approximate  4-year  period  in 
spring  temperatures  may  be  found  at  many  stations  across  the 
province  although  the  true  existence  of  these  periods  at  the 
stations  with  the  low  confidence  level  is  doubtful. 

Closer  scrutiny  of  the  individual  months  used  to  derive 
the  winter  temperature  values  showed  that,  for  December 
temperatures,  all  stations  except  RAN  nad  significant  peaks 
at  the  95%  level  for  a  3.12-  to  3.2  year  period;  that,  for 
January  temperatures,  only  RAN  and  YXD  had  such  peaks  for  a 
3.28-year  period;  and  that,  for  February  temperatures,  six 
stations  had  such  peaks  for  a  3.28-  to  3.39-year  period 
while  the  other  two  stations  had  peaks  nor  a  3.j>6—  to  3. 66- 
year  period. 

A  similar  comparison  done  for  spring  temperatures 
showed  that,  for  March  temperatures,  only  I'VR  and  HAN  die 
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not  have  significant  peaks  for  a  4-year  period;  that,  for 
April  temperatures,  no  station  had  any  significant  peak 
close  to  a  4~year  period;  and  that,  for  May  temperatures, 
only  RAN  and  IXD  had  significant  peaks  for  a  4-year  period. 

7 • 4  Northern  Alberta 

Figure  A. 2  in  Appendix  A  shows  the  location  of  all 
Alberta  stations  used  in  this  study.  BEA  and  FVR  results  are 
discussed  separately  in  this  section.  All  other  stations  are 
dicussed  in  tne  next  two  sections.  The  terms  Northern, 
Central  and  Southern  Alberta  have  been  losely  applied  and 
are  used  here  as  a  convenient  way  of  grouping  the  stations. 

Tables  7.6  to  7.9  present  the  results  obtained  for  FVR 
and  BEA  based  on  the  seasonal,  annual  and  individual  monthly 
values  of  both  temperature  and  precipitation.  Only 
significant  spectral  peaks  (>95%)  are  listed  in  these 
tables.  All  peaks  were  obtained  from  the  unsmoothed 
periodogram  analyses. 

A  marked  feature  found  in  the  analyses  of  FVR  values 
was  the  relative  absence  of  spectral  lines  with  values  less 
than  the  5%  significance  level.  This  result  was  common  to 
most  of  the  stations  analyzed.  The  stations  with  the 
greatest  number  of  such  gaps  analyzed  were  BEA,  YXD  and  YXH . 
This  finding  cast  considerable  doubt  on  the  usefulness  of 
significant  gaps.  Alternatively,  the  question  arose  as  to 
the  meaning  of  a  significant  gap  for  the  purposes  of  this 
study,  since  the  amplitude  of  the  wave  associated  with  such 


Table  7.6.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  annual  and 
seasonal  values  at  FVR. 


P(%) 

Temperatur 

e 

Precipitation 

99 

97.5 

95 

99  97.5  95 

Annual 

128. 

64 . 0 

14.  2 

128. 

16.0 

5.33 

Winter 

16.  0 

128.  42.7 

64.0++  4.27 

Spring 

128. 

16.  0 

6.74 

14.2 

Summer 

32.0 

128. 

8.00 

3.39 

42.7 

2.67 

Fall 

5.  12 

4.92 

none 

indicates 

that  the  peak 

exceeded  the  99.9%  level. 

a  gap  would  be  near  zero.  It  is  of  interest,  for  some 
purposes,  to  know  which  periods  are  not  present  in  a  time 
series,  but,  this  is  not  the  case  for  meteorological  time- 
series  prediction.  For  these  reasons,  the  significant  gaps 
analyzed  were  not  included  in  the  tables  presented  in  this 
chapter . 

Table  7.6  presents  the  analyses  of  annual  and  seasonal 
values  of  temperature  and  precipitation  at  FVR.  Periods  of 
64  and  128  years  appear  a  number  of  times  in  this  table  as 
well  as  in  Tanle  7.7.  Such  periods  have  little  meaning  when 
analyzed  by  the  per iodogram  method.  ilEM  was  used  to  obtain 
an  estimate  of  the  possible  long  period  indicated  by  the 
periodogram  analysis  of  FVR  annual  temperature  values.  The 
analysis  was  performed  using  22  filter  coefficients  on  66 
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Table  7.7.  Periods  (years)  of  significant  peaks  and 
signif icance  levels  (P,%)  for  monthly  values 


at  FVR 

■  9 

i’emperat  are 

Precipitation 

P(%)  99 

97.5  95 

99  97.5 

95 

Jan . 

18.3  16.0 

2.25 

3.39 

Feb. 

3.66 

64.0*  128. 

Mar. 

16.0  14.2 

14.2 

6.40 

Apr . 

3.05 

none 

May 

128.  42.7 

32.0 

5.  82 

Jun . 

32.0  128. 

none 

Jul . 

42.7  32.0 

18.3  5,57 

3.39 

2.29 

Aug . 

32.  0 
8.0  0 

12.  8 

32.0 

14.2 

Sep . 

none 

2.13 

Oct . 

18.3 

none 

Nov . 

2 . 56 

2.61 

128. 

64.0 

5.57 

Pec. 

3.12 

64. 0+  +  128. 

42.7 

5.  12 

+ 

indicates 

indicates 

that  the  peak 
that  the  peak 

exceeded  the  99.5%  level 
exceeded  the  99.9%  revel 

f 

e 

data  values.  MEM  revealed  a  peak  at  the  82.5-year  period. 
However,  the  number  of  data  values  used  and  the  initial 
phase  may  have  shifted  this  peak  from  the  true  value,  as 
described  in  Section  7.1.  An  estimate  of  the  true  period  of 
this  peak  might  be  82. 5±10  years. 


based  on  the  results 
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Table  7.8.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  annual  and 
seasonal  values  at  BEA. 


iemperat  ure 

Precipitat ion 

P  (SJ 

99 

97.5 

95 

99  97.5  95 

Annual 

16.0  + 

64.0 

7.11 

Winter 

16.0 

3.20 

2.21 

none 

Spring 

16.0  + 

none 

Summer 

2.37 

16.  0 
2.  13 

7.  11 

Fall 

4.92 

6.40 

2.56 

32.0  3.05 

+  indicates 

that  the 

peak 

exceeded  the  99.5%  level. 

obtained  during  the  analysis  of  the  test  data.  Such  an 
estimate  was  of  little  use  in  finding  true  periodicities. 

Assuming,  for  the  moment,  that  an  82.5~year  period  did 
exist  in  the  FVR  annual  temperature  varues,  how  real  would 
this  period  re?  Chapter  2  showed  that  there  were  no  monthly 
mean  temperature  values  missing  in  the  EVR  record.  However, 
a  change  of  location  of  the  observing  site  was  noted  in 
Chapter  3.  This  change  took  place  in  1936,  a  few  years  prior 
to  the  mid- point  of  the  data  record,  and  involved  a  move  of 
nearly  8km  to  the  new  site.  This  move  could  have  easily 
resulted  in  a  change  in  the  surrounding  temperature  regime 
which  could  be  interpreted  as  a  portion  of  a  very  long 
period.  The  existence  of  physical  difrerences  between  the 
two  sites  is  not  known  and  could  only  be  determined  by  an 
on-site  investigation. 
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Table  7.9.  Periods  lyears)  of  significant  peaks  and 
significance  levels  (P,%)  for  monthly  values 
at  nEA. 


Temperature 

Precipitat ion 

P {%)  99 

97.5 

95 

99  97.5 

95 

Jan . 

16.0 

none 

Feb . 

3.56 

3 . 56 
3.05 

War .  16.0 

10.7 

4.00 

4.00 

64.0 

16.0 

Apr . 

3.  05 

none 

Way 

64. 0 

none 

Jun . 

2.13 

3.  76 

Jul .  16.0 

2.06 

Aug . 

2.37 

5 .  82 

Sep. 

6.40 

4.92 

2.46 

Oct . 

5.33 

3.05 

2.00 

21.3 

5.33 

Nov.  2.56++ 

32  «  0  + 

6.40 

Dec . 

3.20 

2.21 

7.11 

8.00 

+  indicates 
+  +  indicates 

that  the 
that  the 

peak 

peak 

exceeded  the  99.5%  level, 
exceeded  the  99.9%  level. 

An  article 

recen  tly 

written  by  Schaal  and  Dale 

(1977) 

creates  more  doubt  about  the  existence  of  long  periods  in 
data  samples.  Briefly,  they  found  that,  because  of  a  change 
in  the  time  period  used  to  describe  the  "climatological 
day",  the  mean  March  temperature  over  Indiana  acquired  a 
bias  of  0.67C  which  is  nearly  equal  to  the  amount  or 
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climatic  ,fcooiingu  observed  in  the  mean  March  temperature 
values  of  Indiana  over  the  last  40  years.  Mr.  T.  Donnelly, 
Chief  of  Surface  Inspectors  at  WRO,  related  that,  shortly 
after  the  second  World  War,  the  World  Meteorological 
Organization  standardized  the  observation  times.  The  end  of 
the  climatological  day  in  Canada  was  moved  from  0000GMT  to 
0600GMT  for  tne  AES-controled  stations.  Mr.  Donnelly  also 
advised  that  no  regulations  are  imposed  at  privately- 
operated  stations  and  that  the  time  of  observation  may  vary, 
presumably  as  the  work.  load  of  the  observer  varies. 

In  view  of  these  findings,  it  is  impossible  to  say 
whether  a  period  of  82.5  years  with  an  amplitude  of  0.64C  at 
FVR  is  real  or  results  from  some  or  all  of  the  factors 
mentioned  above.  This  uncertainty  added  to  the  reluctance  to 
use  MEM  because  apparent  waves  could  have  been  introduced 
into  the  data  by  one  or  more  obscure  factors. 

It  is  interesting  to  note  that  the  results  for  BEA, 
shown  in  Tables  7.8  and  7.9,  contained  few  very  long 
periods.  The  site  at  BEA  was  moved  once  in  1958. 
Simultaneous  records  were  kept  at  the  two  sites  during  the 
three-year  period  preceding  the  move,  carder  (1962)  referred 
to  this  move  and  stated  that  the  comparison,  which  he 
carried  out,  aid  not  reveal  an y  significant  difference  in. 
the  temperature  and  precipitation  records.  Unfortunately, 
the  data  collected  during  these  three  years  were  never 
published  and  are  not  available  from  the  station.  However , 
accepting  Carder* s  statement  would  be  eguivalent  to  ignoring 
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the  change  of  location  of  the  observing  site, 
result  affected  by  such  an  assumption  would  be  rhe  64-year 
period  found  an  the  annual  temperature  values.  However,  too 
many  unknown  variables  remain  and  any  conclusion  associated 
with  this  64-year  period  could  only  be  called  irresponsible. 

All  further  discussion  in  this  chapter  deals  with 
significant  peaks  with  periods  of  less  than  32  years. 

Very  few  analyses  of  FVR  or  BEA  data  values  showed  any 
"unexpected"  results  because,  as  mentioned  in  Section  7.2, 
5%  of  the  analyzed  spectral  lines  were  expected  to  exceed 
the  95%  significance  level. 

The  fairly  frequent  occurrence  of  a  16-year  peak, 
mainly  in  the  temperature  analyses  for  FVR  and  BEA,  was 
investigated  but  little  agreement  was  round  among  the  phase 
angles.  Regardless  of  the  99.5%  significance  level  achieved 
by  the  16-year  period  in  the  analyses  of  annual  and  spring 
temperature  values  at  BEA  (see  Table  7.8),  there  was 
insufficient  evidence  available  for  any  positive  statement. 

The  shorter  periods,  towards  the  2-year  "Nyguist" 
period,  showed  little  organization.  Too  few  of  these 
periods,  derived  from  the  analysis  of  individual  months, 
were  evident  in  the  results  of  analyses  of  seasonal  or 
annual  data.  There  was  no  pattern  evident,  either  at  FVR  or 
BEA,  or  common  to  both,  other  than  that  discussed  in  the 
previous  section. 

The  data  values  which  were  listed  as  missing  in  Table 
2.2,  were  or  little  concern  because  few  of  the  results  of 
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the  precipitation  analyses  were  of  interest.  The  two  values, 
missing  from  the  BEA  record,  were  for  August  1915  and 
January  1916 ;  both  months  were  in  the  first  year  of 
observations.  The  five  missing  values  for  the  FVE 
precipitation  record  were  September-Cctober  1911,  March- 
Aprii  1916  and  April  1934. 

7.5  Central  Alberta 

Tables  7.10  to  7.15  present  the  results  obtained  for 
HAN,  YXD  and  LAC  based  on  the  seasonal,  annual  and 
individual  monthly  values  of  both  temperature  and 
precipitation.  Only  significant  spectral  peaks  (>95%) , 
obtained  from  the  unsmoothed  periodogram  analyses,  are 
listed  in  these  tables. 

The  three  stations  considered  in  this  section  exhibit 
patterns  similar  to  those  which  have  axready  been  discussed 
in  Section  7.3.  Periods  of  the  order  of  2.57  years  appear  in 
a  number  of  places  but  with  insufficient  regularity  and 
significance  to  indicate  extraordinary  circumstances.  There 
are,  however,  one  or  two  periods  worthy  of  mention. 

First,  consider  the  18.3“  to  21.3-year  period  analyzed 
from  the  July  temperatures  of  both  BAN  and  LAC.  The  18.3- 
year  period  was  also  apparent  in  the  January  temperatures  of 
these  two  stations.  Of  the  four  occurrences  of  this  peak, 
only  the  July  peaks  for  RAN  exceeded  the  expected  number  of 
three  significant  spectral  peaks.  The  January  peak  for  LAC 
lacked  credibility  because  it  was  the  only  spectral  line 
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Table  7.10.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  annual  and 
seasonal  values  at  RAN, 


Temperat  ure 

Precipitation 

99  97.5 

95 

99 

9  7 . 5 

95 

Ann  ual 

none 

none 

Winter 

3.28 

2.25 

2.37 

Spring 

3.88 

2.42 

128. 

64.0 

Summer 

42.7 

21.3 

none 

Fall 

5.  12 

4.92 

2.84 

2.78 

which 

exceeded  the 

95% 

signified  nee 

level . 

A  few 

observations  are  presented  concerning  the  July  peaks. 

The  weignted  frequencies  indicated  significant  periods 
of  20.04  and  19.78  years  for  RAN  and  LAC,  respectively.  The 
phase  angle  for  the  RAN  peak  was  found  to  be  41°,  based  on 
1905,  with  an  amplitude  of  0.66C.  Tne  phase  angle  for  the 
LAC  peak  was  found  to  be  104°,  based  on  1908,  with  an 
amplitude  of  0.64C;  the  phase  angie,  when  adjusted  by 
1 8°/year  (based  on , an  exact  20-year  period},  became  50°. 
However,  these  results  are  not  necessarily  as  impressive  as 
they  seem  because  YXD,  which  is  approximately  equidistant 
from  RAN  and  LAC,  showed  no  evidence  of  this  period  in  the 

July  temperatures. 

The  movements  of  the  observing  sites  at  LAC  and  RAN 
offered  little  help  in  identifying  a  reason  which  would  have 
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Table  7.11.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  monthly  values 
at  RAN. 


Temperature  Precipitation 


P  {%)  99 

97.5 

95 

99  97.5 

95 

Jan . 

18.3 

3.28 

2.25 

none 

Feb. 

3.28 

18,3 

21.3 

2.  10 

Mar . 

3.88 

4.00 

8.00 

2.46 

Apr. 

none 

11.6 

2.61 

May 

42.  7 
4.00 

none 

Jun . 

12.  8 

2.37 

Jul . 

64.0 

21.3 

9.  85 

18.  3 
9.85 

3 . 56 

Aug. 

none 

none 

Sep . 

4.92 

none 

Oct . 

none 

2.67 

5.33 

2.84 

Nov.  2.56 

2.61 

5.12 

2.25 

42.7 

Dec. 

/ 

none 

2.  25 

8.00 

2.21 

explained  thus 

period 

,  especially  since  the  LAC  site 

was  in 

the  same  place 

for  65 

of  the  68 

data  years  used.  The 

missing 

data  offered 

little 

help  as 

well  because  the 

only 

temperature  value  missing 
for  November  1913  at  LAC. 


from  the  RAN  and  LAC  data  sets  was 
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Table  7.12.  Periods  (years)  of  significant  peaks  and 
significance  levels  (?,%)  for  annual  and 
seasonal  values  at  YXD. 


Temperature 

Precipitation 

93  97.5 

95 

99 

97.5  95 

Annual 

none 

10.7  2.25 

Winter 

3.28 

10.  7 
2.21 

3.20 

Spring 

5.12 

4.00 

2.61 

none 

Summer 

128.*-+ 

32.  0 

10.7 

12.8 

Fall 

4.92  + 

4.  57 

+  indicates  that  the 
++  indicates  that  the 

peak 

peak 

exceeded 

exceeded 

t  he 
t  he 

99.5%  level, 

99=9%  level. 

Very 

little  infor 

mation 

could 

be 

obtained  from  the 

analyses  of  precipitation  values  at  the  Central  Alberta 
stations.  The  most  noteworthy  result  was  a  10. 7-year  period 
obtained  from  the  analysis  of  summer  precipitation  values 
for  YXD.  Thrs  period  was  also  the  most  significant  peak  in 
the  YXD  annual  precipitation  analysis.  This  result  could 
hardly  be  considered  as  overwhelming.  However,  a  brief  test 
was  performed  on  the  summer  peak  because  this  peax.  accounted 
for  9.2%  of  tne  variance  in  the  summer  data  while  the  same 
peak  in  the  annual  values  accounted  for  8.5%  of  the 

variance. 

Data  were  generated  for  a  single  10.7~year  wave  with  a 
mean  of  227.4mm,  amplitude  of  28.7mm  and  phase  angle  of 
77.66°  based  on  an  origin  of  1882.  These  values  were 
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TaDle  7.13.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  monthly  values 
at  IXD. 


P(%) 

Temperat  ure 

Precipitation 

99  97.5 

95 

99  97.5 

95 

Jan . 

2.  13 

3.28 

none 

Pen. 

3.28 

4,27 

2.33 

Mar . 

4.00 

10.7 

2.33 

4.00 

Apr . 

3.05 

2.72 

3.  66 

6.74 

2.61 

May 

4.00 

11.6 

8.53 

3.88 

4.  13 

Jun . 

2.91 

128. 

9.85 

2.98 

12.8 

2.  13 

6.40 

4.74 

Jul . 

128** 

2.03 

5.33 

11.6 

Aug . 

128. 

4.27 

10.7 

3.  88 

3.  12 

Sep . 

9.85 

1 4 . 2 

6.40 

11.6 

4.92 

5.33 

4.92 

Oct . 

6.40 

2. 03  2.67 

2.91 

Nov . 

2.56  + 

4.57 

32.0 

4.92 

2.61 

Dec . 

3.  12 

2.21 

2.25 

2.46 

2.21 

+  indicates  that  the  peak  exceeded  the  99.5%  level, 
+  -<-  indicates  that  the  peak  exceeded  the  99.9%  level. 
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Table  7,14.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  annual  and 
seasonal  values  at  LAC. 


Temperat  ure 

Precipitation 

p  m 

99  97 .5  95 

99  97.5  95 

Annual 

2.67 

2.61 

9.85 

9.  14 

Winter 

none 

4. 74 

Spring 

3.88 

none 

Summer 

25.  6 

3.80 

fall 

5.  12 
4.92 

10.  7 

6.40 
2. 51 

obtained  froo  the  analysis  of  summer  precipitation  data.  The 
sums  of  the  squares  of  the  differences  between  the  actual 
and  generated  values  and  between  the  actual  values  and  the 
data  mean  were  calculated,  i.e.,  using  a  persistence 
forecast  in  tne  second  case.  The  sum  obtained  using  the 
generated  values  was  13%  smaller  than  the  sum  obtained  using 
persistence.  This  would  seem  to  indicate  some  measure  of 
validity  in  assuming  the  presence  of  a  10.7-year  period. 

The  smoothed  per iodogram  provided  slightly  different 
results  than  the  raw  periodogram.  The  raw  periodogram 
derived  10.7-  and  12.8-year  peaks  from  the  YXD  summer 
precipitaion  data.  These  peaks  were  separated  by  a  non¬ 
significant  11.6-year  period.  The  smoothed  periodogram 
(using  a  three-point  Daniell  filter)  derived  five 
significant  peaks  among  which  the  11. 6- year  peak  exceeded 
the  99.5%  significance  level.  Difficulties  arose  in  trying 
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Table  7.15.  Periods  lyears)  of  significant  peaks  and 
significance  levels  (P,%)  for  monthly  values 
at  LAC. 


Temperature 

Precipit at  ion 

P(») 

99  97.5 

95 

99  97.5 

95 

Jan . 

18.3 

2.84 

Peb . 

3.28 

3.05 

Har . 

3.88 

3.88 

Apr. 

3.05 

2.61 

2.56 

Hay 

none 

9.14 

Jun . 

none 

14.2 

Jul . 

64.0  21.3 

18.3 

3.66 

Aug. 

8.00 

9.85 

5.  12 

Sep . 

4.92 

none 

Oct . 

none 

9. 14 

Nov . 

2.56*  2.61 

4.92 

2.  56 

Pec . 

3.12 

none 

*  indicates  that  the 

peak 

exceeded  the  99.5%  level 

e 

to  interpret  the  results  of  the  smoothed  pericdogram  because 
these  five  peaks  were  no  longer  independent.  The  number  of 
degrees  of  freedom  used  to  calculate  the  confidence  interval 
was  approximately  4.5  for  the  smoothed  periodogram  as 
opposed  to  1.5  for  the  unsmoothed  periodogram.  However, 
unlike  the  ACV  method  where  a  smalr  lag  value  yields  few 
spectral  lines,  the  modified  periodogram  analyzes  the  same 
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number  of  spectral  lines  for  both  smoothed  and  unsmoothed 
results.  How  Lhen  are  the  five  adjacent  peaks  assessed  and 
what  are  the  implications  of  smoothing?  These  were  questions 
which  would  nave  proven  very  difficult  to  answer  and,  based 
on  the  results  obtained  while  using  the  smoothed 
periodogram,  would  likely  have  resulted  in  little  additional 
information  than  was  available  from  the  unsmoothed 
periodogram. 

7*6  Southern  Alberta 

Tables  7.16  to  7.21  present  the  results  obtained  for 
1TH,  YXH  and  HAN  based  on  the  seasonal,  annual  and 
individual  monthly  values  of  both  temperature  and 
precipitation.  Only  significant  spectral  peaks  i>95%) , 
obtained  from  the  unsmoothed  periodogram  analyses,  are 
listed  in  these  tables. 

Neglecting  the  periods  already  discussed  in  Section 
7.3,  there  remain  one  or  two  points  worthy  of  mention 
concerning  these  last  three  stations. 

YXH  was  the  only  station  which  exhibited  an 
"unexpected"  number  of  significant  peaxs  on  more  than  one 
occasion.  Some  of  the  occurrences  invoxved  very  long  periods 
which  were  readily  dismissed  owing  to  the  fact  that  the 
station  was  moved  out  of  the  river  valley  in  1930.  However, 
there  renrained  the  frequent  occurrence  of  an  18.3-  to  21.3- 
year  period  in  the  precipitation  data,  as  shown  in  Tables 


7.18  and  7.19. 
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Table  7.1t>.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  annual  and 
seasonal  values  at  LTH . 


Temperat  ure 

Precipitation 

p(%) 

99  97.5  95 

99  97.5  95 

Annual 

2.67 

12.8 

4.27 

Winter 

6.  74 

3.  28 

64.0 

Spring 

3.88 

none 

Summer 

12.8 

42.7 

12.  8 

Fall 

5.  12 
4.92 

7.53 

none 

A  true  wave  was  generated  using  a  mean  of  an 


amplitude  of  17.3mm,  a 

period 

of  19.7  years. 

and 

a 

phase 

angle  of  114°  based  on 

an  origin  in  1884.  The 

am  pi 

itude  was 

obtained  by  taking  the 

sum  of 

the  amplitudes 

of 

the 

two 

peaks  present  in  the  analysis  of  winter  precipitation  at 
YXH ;  this  procedure  was  based  on  the  assumption  that  the  two 
peaks  were  located  on  either  side  of  the  true  significant 
peak.  Two  sums  were  calculated  in  a  manner  analogous  to  the 
calculations  performed  on  the  YXD  data  in  the  previous 
section.  The  rirst  sum,  obtained  by  using  the  generated  wave 
and  the  winter  precipitation  data,  was  equal  to  73%  of  the 
second  sura  obtained  using  persistence.  Again,  the  analyzed 
peak  seemed  to  indicate  some  measure  of  fit  when  applied  to 
the  data.  However,  it  should  be  noted  that  persistence  was 
much  better  in  the  last  30  to  35  years  of  the  data  record. 
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Table 

7.17.  Periods 
significance 
at  LTH. 

(years)  of 
levels  (P/%) 

significant  peaks  and 
for  monthly  values 

Temperat  ur 

e 

Precipitation 

p  i 

99  97.5 

95 

99  97.5  95 

Jan . 

6.74 

64.0 

7.11 

Feb . 

3.  28 

9.85 

32.0 

2.21 

Mar . 

3.  88 

2.78 

2.25 

2.21 

Apr . 

2.78 

2.61 

May 

2.00 

3.66 

Jun . 

12.8 

none 

Jul . 

3.39 

3.56 

Aug . 

12.  8 

2.  72 

8.00 

Sep . 

7.53 

none 

4.92 

Oct . 

none 

none 

Nov . 

2.56 

2.61 

2.06 

Dec. 

3.  12 

4.4  1 

2.46 

The  reason  for  this  may  possibly  be  seen  by  considering  the 
data  plot  in  figure  7.8.  The  unsmoothed  periodogram  analysis 
of  these  data  was  shown  in  Figure  7.3. 

The  20-year  period  can  be  seen,  to  a  certain  extent, 
between  1884  and  1930.  Presumably,  this  early  existence  of 
the  wave  was  sufficient  in  order  to  be  analyzed  by  the 
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Table  7.18.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P,%)  for  annual  and 
seasonal  values  at  YkH. 


Temperature 

Precipitation 

P  (%)  99 

97.5 

95 

99  97.5 

95 

Annual  64.0* 

128. 

4.9  2 

18.  3 

25.6 

14.2 

2.37 

Winter 

3.28 

2.46 

21.3*  2.78 

128. 

18.3* 

Spring 

128. 

25.6 

64.  0 

18.3 

4.00 

6.74 

2.61 

Summer  64.0** 

128. 

2.37* 

12.8 

42.  7 

12.8 

4.00 

Fall  4.92* 

7.53 

16.0 

*  indicates 

that  the 

peak 

exceeded  the  99.5%  level 

r 

**  indicates 

that  the 

peak 

exceeded  the  99.9%  level 

• 

periodograra  method.  It  is  conceivable  that  this  wave  was  the 
result  of  tne  numerous  changes  of  location  of  the 
meteorological  site  rn  the  period  between  1883  and  1930.  The 
station  remained  approximately  10  years  at  each  of  the  first 
three  sites  and  20  years  at  the  fourth  site.  The  station  has 
not  been  moved  since  1938.  Based  on  this  information,  it 
becomes  impossible  to  state  whether  the  20-year  period  is 
truly  a  climatic  fluctuation  or  results  from  the  many 
relocations  or  the  station. 

Slutzky  (1937),  using  a  rigorous  mathematical  approach, 
demonstrated  that  series  of  random  data  could  provide 


Table  7.19.  Periods  (years)  of  significant  peaks  and 
significance  levels  ( Pt% )  for  monthly  values 
at  YXH. 


Temperat  ure 

Precipitation 

P(%)  99  97.5 

95 

99  97.5 

95 

Jan . 

2.13 

21.3 

128. 

18.3 

3.05 

2.78 

Feb. 

3.28 

21.3**  42.7 

18.3 

2.46 

War . 

4.00 

18.3*  64.0 

42.7 

2.17 

14.2 

Apr. 

128. 

6.74 

128. 

3.05 

3.56 

Way 

2.  17 

3.56 

42.7 

Jun . 

64.0*  42.7 

2.  37* 

2.  8a 

12.8 

Jul . 

64.0*  3.39 

128. 

5.  57 

2.03 

3.88 

Aug . 

12.8 

8.00 

4.00 

6.74 

25.6 

Sep. 

4.92 

8.00 

Oct . 

2.91 

18.3 

18.3 

5.33 

Nov . 

4.92 

2.56 

16.0 

4.92 

Dec. 

2.  46 

3.12 

8.00 

3.56 

2.  78 

indicates  -chat  the 

peak 

exceeded  the  99.573  level 

t 

+t 

indicates  that  the 

peak 

exceeded  the  99.3%  level 

• 

evidence  of  periodic  fluctuations  either  throughout  the 
whole  series  or  in  segments  of  the  series.  These  findings 
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Table  7.20.  Periods  (years)  of  significant  peaks  and 
significance  levels  (P ,%)  for  annual  and 
seasonal  values  at  MAN. 


P(%) 

Temperat  ur e 

Precipitation 

99  97.5 

95 

99  97.5  95 

Annual 

21.3 

64.0 

Winter 

3.20 

16.0 

64 . 0  + 

7.11 

3.39 

Spring 

2.  46 

10.7 

Summer 

none 

none 

Fall 

4.92 

4.57 

9.14 

+  indicates  that  the 

peak 

exceeded  the  99.5%  level. 

YEAR 

Figure  7.8.  YXH  winter  precipitation. 

were  especially  valid  for  data  series  which  did  not  have 
totally  independent  elements.  The  analogy  of  this  last 
finding  fo  meteorological  data  is  obvious,  even  if  the  data 
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Table 

7.21.  Periods 
significance  1 
at  MAN. 

(years)  of 
evels  (P,%) 

significant  peaks  and 
for  monthly  values 

Temper at ure 

Precipit  at ion 

P(%) 

99  97.5 

95 

99  97.5  95 

Jan . 

7.1 1 

64. 0 

feb . 

10.7  3.39 

none 

Mar . 

none 

12.8 

2.56 

Apr . 

none 

none 

May 

none 

10.7  3.20 

Jun . 

12.  8 

2.91 

Jul . 

64.0 

3.39 

3.39 

Aug . 

10.  7 

8.00 

8.00 

Sep. 

none 

3.  76 
3.39 

Oct . 

21.3 

none 

Nov . 

2.56 

4.92 

9.14 

Dec . 

16.0 

3.20 

none 

consist  of  mean  temperatures  for  a  particular  month  or 
season. 

The  changes  of  the  sites  at  LTH  and  MAN  were  not  as 
numerous  nor  as  drastic  as  the  changes  at  YXH.  Similarly, 
the  results  of  the  analyses  for  LTH  and  MAN  were  not  as 
dramatic  as  tnose  for  YXH. 

Missing  data  in  the  YXH  record  may  have  contributed 
also  to  the  generation  of  erroneous  results.  The  winter 
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precipitation  data  contained  three  values  which  had  been 
estimated,  wnolly  or  in  part,  or  obtained  from  the  original 
published  records.  However,  the  information  in  these 
original  records  was  rejected  in  compiling  the  microfiche 
data  files  which  were  described  in  Chapter  2.  Two  of  the 
three  estimated  values  were  varied  by  50mm.  The  resultant 
spectrum  revealed  the  same  significant  peaks  as  previously 
mentioned.  This  would  seem  to  indicate  that  the  existence  of 
three  estimated  values  in  this  data  set  did  not  contribute 
significantly  to  the  resultant  spectrum.  The  missing  values 
of  monthly  precipitation  were  October  1883,  November- 
December  1885,  January-April  1886,  January -March  1887, 
August  1889,  October  1889,  November-December  1910  and 
January  1911. 

No  other  periods  were  felt  to  be  truly  significant. 


CHAPTER  8 


SUMMARY  AND  CONCLUSIONS 

8 . 1  Spectral  Analysis  Methods 

The  set  of  programs  developed  for  this  study  provided 
the  capability  of  using  the  modified  periodogram  method  and 
the  Maximum  Entropy  Method.  The  periodogram  could  calculate 
raw  and  smoothed  values  of  phase,  frequency  and  amplitude 
squared.  Smoothing  was  performed  by  a  Daniell  filter. 
Tapering  of  the  ends  of  the  data  was  available.  Three 
subroutines  provided  an  interface  between  the  spectral 
output  and  a  slightly  modified  version  of  a  plotting 
subroutine  available  from  the  Computing  Services  Program 
Library  of  the  University  of  Alberta. 

The  unsmoothed  modified  periodogram  method  was  used  in 
all  cases  studied.  It  was  found  to  be  a  highly  efficient 
analysis  tecnnique  and  it  performed  excellently  on  the  test 
data  which  were  different  combinations  of  cosine  waves. 

The  smoothed  phase  anqle  and  the  weighted  frequency 
proved  to  be  quite  useful  in  providing  better  estimates  of 
the  phase  angxe  and  frequency  than  the  unsmoothed  values. 
Tapering  was  round  to  offer  no  advantages.  Caution  should  be 
exercised  in  adopting  these  features  of  the  program. 
Analyses  of  test  data  would  provide  a  better  understanding 
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of  its  capabilities. 

The  smoothed  modified  periodogram  was  calculated  for 
all  data  groupings.  Smoothing  did  add  a  measure  of  stability 
to  the  results,  but,  for  the  purposes  of  this  study, 
provided  no  improvement  on  the  information  available  from 
the  unsmoothed  periodogram. 

Robinson  (1967}  describes  how  the  Subroutine  NLOGN  (see 
Section  4.3  and  Table  B.8)  could  be  used  to  calculate  two 
real-valued  signals  simultaneously,  and  to  calculate  the 
cross  correlation  between  two  real-valued  signals.  Both  of 
these  modifications  would  be  desirable  additions  to  the 
program  prepared  for  this  study. 

The  Maximum  Entropy  Method  was  used  very  little.  This 
was  due  in  part  to  the  problems  of  interpreting  the  results, 
as  described  in  Section  7.1.  Also,  there  are  the 
difficulties  presented  by  inhomogeneous  data.  Results 
obtained  by  using  MEM  acquire  dubious  credibility  given  the 
discussion  on  the  existence  of  long  waves,  presented  in 
Chapter  7.  This  would  be  true  even  if  MEM  were  improved 
beyond  the  technique  presented  by  Ulrych  and  Clayton  (1976). 
The  greater  improvement  required  is  in  the  reliability  of 
the  meteorological  data. 

8.2  Climatological  Analysis 

The  data  used  in  this  study  were  considered  to  be  as 
good  as  any  available.  However,  findings  such  as  the  ones 
presented  by  Schaal  and  Dale  (1977)  illustrate  the  degree  of 
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inhomogeneity  of  the  data. 

Comprehensive  histories  of  the  meteorological  stations 
are  sadly  lacking.  The  compilation  of  such  histories  for 
only  eight  stations  demanded  much  more  time  than  should  have 
been  necessary,  but  these  histories  were  required  in  order 
to  evaluate  tne  results  properly.  The  histories  presented  in 
Chapter  3  do  not  contain  all  the  information  which  should  be 
considered  when  making  conclusions  but,  as  indicated  in 
Chapter  1 ,  changes  of  instrumentation,  observer  or  observing 
practices  are  frequently  difficult  if  not  impossible  to 
ascertain. 

Tiiree  observations  may  be  made  following  the 
climatological  inves tigations .  These  are  the  three  distinct 
periods  revealed  by  the  analysis  of  seasonal  temperatures. 

The  most  notable  of  the  three  observations  involves  the 
five-year  period  revealed  in  the  fall  temperatures  at  all 
stations.  No  explanation  can  be  formulated  at  this  time,  to 
dispel  the  existence  of  this  period.  Little  evidence  was 
found  to  indicate  that  the  period  was  the  result  of  leakage 
of  power  from  low  frequencies  or  aliasing  of  high 
frequencies.  The  occurrence  of  this  period  at  all  stations 
and  the  similarity  of  the  phase  angles  supported  the  theory 
that  this  period  may  be  guite  real. 

Georgiades  (1977)  resolved  a  five-year  period  for 
annual  mean  temperatures  at  Regina  and  Saskatoon, 
Saskatchewan.  Only  the  annual  mean  temperatures  for  YXH  gave 
evidence  of  a  five-year  period  (see  Table  7. 1) •  No  reasons 
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can  be  presented  at  this  time,  which  would  explain  this 
apparent  difference  between  Saskatchewan  and  Alberta. 

The  two  remaining  observations  are  not  as  notable  as 
the  first.  However,  there  appears  to  be  evidence  of  fairly 
widespread  tnree-  and  four-year  periods  in  the  winter  and 
spring  mean  temperatures,  respectively. 

Closer  investigation  of  these  tnree  periods  should  be 
made  in  order  to  determine  their  possible  causes  and  the 
extent  of  their  existence  across  Canada. 
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APPENDIX  A 


TOPOGRAPHICAL  HAPS 

This  Appendix  contains  a  set  of  maps,  the  first  of 
which  shows  tne  geographical  distribution  across  Alberta  of 
the  eight  stations  used  in  this  study  (Figure  A. 2) ;  all 
other  maps  show  the  detailed  geography  at  each  location 
(Figures  A .3  to  A. 10). 

The  detailed  maps  were  constructed  using  the 
topographical  series  available  from  the  Surveys  and  Happing 
Branch,  Department  of  Energy,  Hines  and  Resources.  All 
elevations  are  given  in  feet  and  can  be  convened  using  the 
scale  shown  in  Figure  A.1. 
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Figure  A.1.  Conversion  scale  for  elevations. 
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Scale  1:50,000  contour  interval  so  feet 

All  Elevations  in  Feet  above  Mean  Sea  Level 
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2  -  (1959-1975) 


Figure  A. 3.  Beaverlodge  CDA  site  locations 
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Edmonton  City  Site  locations. 
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Figure  A- 5.  Fort  Vermilion  C DA  site  locations 
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Scale  1 :50,000  Echelle 
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Figure  A.  6.  Laconibe  CDA  site  locations. 
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Figure  A. 7.  Lethbridge  CDA  site  locations. 
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1  -  (1928-1944) 

2  -  (1944-1949) 

3  -  (1949-1975) 


Figure  A. 8.  Manyberries  CDA  site  locations. 
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Figure  A. 9.  Medicine  Hat  site  locations. 


I 
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Figure  A. 10.  Ranfurly  site  iocations. 


APPENDIX  B 


COMPUTES  PROG BA MS  AND 
SAMPLE  PRINT-OUTS 

This  Appendix  contains  a  copy  of  the  computer  programs 
used  in  this  study  (Tables  B.1  to  B.14). 

Sample  print-outs  of  the  FFT  and  MEM  subroutines  are 
also  included  (Tables  B.15  to  B.19).  All  print-outs  were 
obtained  from  the  data  values  listed  in  Table  B.15. 

The  programs  were  written  for  a  FORTRAN  G  compiler. 
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Table  B.l.  Main  program,  as  described  in  Chapter  6. 
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OF  DATA  ON  FIRST  CARD(S)  IF  DIFFERENT  FROM  SUBSEQUENT 


IT  MUST  BE  A  NUMBER 


MAIN  PROGRAM 

***********  *****************#**********i,*£l([;(t;£:$[;<c*  ******************* 

FOLLOWING  ARE  THE  PARAMETERS  TO  BE  INPUT: 

II  =  NUMBER 
CARDS. 

=  0  IF  ALL  CARDS  ARE  THE  SAME. 

NCAS  =  NUMBER  OF  DATA  POINTS,  MAXIMUM  4096. 

NX  =  THE  INTEGER  POWER  OF  2  USED  TO  GET  NZERO ; 

FROM  1  TO  12. 

NZERO  =  TOTAL  NUMBER  OF  DATA  POINTS  NEEDED  TO  HAVE  AN  INTEGER 
POWER  OF  2.  IT  MUST  BE  GREATER  THAN  OR  EQUAL  TO  NCAS;  IT 
MUST  BE  ONE  OF  :  2,4,8,16,32,64,128,256,512,1024,2048,4096. 

JOB  =  JOB  NUMBER  FOR  IDENTIFICATION. 

KP  (POSITIVE)  =  NUMBER  OF  DATA  VALUES  PLOTTED  AND  PRINTED. 

(NEGATIVE)  =  NUMBER  OF  DATA  VALUES  PRINTED  ONLY. 

=  0  NO  DATA  PRINTED  OR  PLOTTED. 

NOTE:  UNLESS  ALTERED,  THE  PLOTTING  SUBROUTINE  PLTG  WILL  REJECT  THE 
JOB  IF  MORE  THAN  1200  DATA  POINTS  ARE  PASSED  ON  TO  IT. 

LB  =  NUMBER  OF  SPACES  ADDED  AT  '’"HE  BEGINNING  OF  THE  DATA 

DISPLACING  THE  FIRST  VALUE  ON  THE  DATA  PLOT  WHICH  ALLOWS  THE 
USE  OF  A  COMMON  ABSCISSA  FOR  PLOTS  WITH  DIFFERENT  FIRST 
ABSCISSA  VALUES.  (NOTE:  NCAS+LB  MUST  BE  LESS  THAN  4097.) 

NOM  =  0  NO  TAPERING  OF  DATA. 

=  NUMBER  OF  DATA  POINTS  TAPERED  AT  EACH  END  (INCLUDING  THE 
FILTER  WEIGHTS  ZERO  AND  ONE),  MAXIMUM  100. 

KT  =  +1  OR  -1  FFT  SPECTRAL  ANALYSIS  ONLY. 

=  + 2  OR  -2  MEM  SPECTRAL  ANALYSIS  ONLY. 

FFT  AND  MEM. 

SPECTRAL  OUTPUT  PLOTTED. 

NO  PLOTTING  OF  SPECTRAL  OUTPUT. 

KN  =  +1  U N3MOOTH ED  FFT  POWER  PLOTTED  ONLY. 

SMOOTHED  FFT  POWER  PLOTTED  ONLY. 


OR 

=  +3  OR  -3 
(POSITIVE) 
(NEGATIVE) 
=  +1 
=  +2 


=  +3  BOTH  -  UNSMOOTHED  AND  SMOOTHED  POWER  PLOTTED. 

(NOTE:  KN  IS  NOT  USED  WHEN  KT  IS  NEGATIVE) 

NOTE:  UNLESS  ALTERED,  THE  PLOTTING  SUBROUTINE  PLTG  WILL  REJECT  THE 
JOB  IF  MORE  THAN  1200  POWER  VALUES  ARE  PASSED  ON  TO  IT. 

NFIL  =  0  NO  FILTERING  OF  F.F.T.  RESULTS. 

=  NUMBER  OF  POINTS  IN  THE  DANIELL  FILTER  WHICH  MUST  BE  PLUS 
OR  MINUS  3,  5,  7,  9  OR  11 . 

(POSITIVE)  FOR  FILTERING  ON  COSINE,  SINE  AND  FREQUENCY. 
(NEGATIVE)  FOR  FILTERING  ON  POWER  SPECTRUM  AS  WELL  AS  THE 
OTHER  THREE  PARAMETERS. 

MEMN  =  NUMBER  OF  FILTER  COEFFICIENTS  FOR  MEM  SPECTRAL  ANALYSIS, 
MAXIMUM  (NCAS-2) . 

=  0  r  NUMBER  OF  COEFFICIENTS  CHOSEN  AUTOMATICALLY  USING 
AKAIKE'S  FINAL  PREDICTION  ERROR. 

MEMX  =  MULTIPLICATION  FACTOR  FOR  EXPANDING  THE  LOWER  END  OF  THE 

FREQUENCY  SPECTRUM.  SET  EQUAL  TO  ONE  IF  NO  EXPANSION  DESIRED. 
SIGN  =  +1.0  OR  -1.0  DETERMINES  THE  SIGN  OF  THE  EXPONENTIAL 
FUNCTION  OF  THE  TRANSFORM  OF  THE  FFT. 

TIME  INTERVAL  OF  THE  DATA  IN  UNITS  CHOSEN  FOR  THE  PROBLEM. 

AND  FOR2  ARE  VARIABLE  FORMATS  (MAX.  60  CHARACTERS) : 

FORI  READS  THE  FIRST  II  VALUES  (OMIT  IF  II  =  0)  , 

FOR2  READS  ALL  OTHER  VALUES. 

******* ************************************** ******* ****** ******** 


DT  = 
FORI 


108 


100  F0RMAT(13I5,F5.1,F5.3) 

101  FORMAT (1  H  1 , ' *****  NX  AND  NZERO  ARE  NOT  COMPATIBLE  ******) 

102  FO  RM  AT  ( 1  5  A4 ) 

103  FORMAT (1H1 , 'JOB  NUMBER',16,*  WITH‘,15,*  VALOES  USED  OF  WHICH1, 15,* 
*  VALOES  ARE  PLOTTED  AND/OR  PRINTED;  THE  SIGN  FOR  THE  FFT  IS',F5.1, 
*///) 

104  FORM AT ( 1  HO  , ' RAW  DATA  MEAN  =  *rFl0.3,'  RAW  DATA  VARIANCE  =  ',F13.3, 
*®  RAW  DATA  STANDARD  DEVIATION  =.',F10.3,/) 

105  FORMAT  (1H0, 1  POINT  NO.  ’,1218) 

106  FORMAT  OH  ,  '  VALUE*  ,  6  X  ,  12F8.2) 

107  FORMAT(lH1,' *****  YOO  BLEW  IT,  KT  =  0  *****') 

Q  ******  *********#***##  jt*#***  ***#*&  ###*## '£*****&#*  ************  ****** 

f 

DIMENSION  DATA  (4  096)  ,FOR1  (15)  ,F0R2(15)  ,V(4096)  , W (  4096)  , PLOT (4 096) 

5  READ (5, 100, END=95 ) II , NCAS , NX , NZ ERO , JOB , K P, LB, NOM , KT , K N, NF IL, MEMN , 
*MEMX , SIGN , DT 
NZO=  2**N  X 

IF  (NZO.EQ. NZERO)  GO  TO  10 
WRITE  (6 , 101) 

STOP 

10  IF  (II. EQ.0) GO  TO  15 
READ  (5, 102)  FORI 
15  READ  (5, 102)  FOR2 

IF  (II.EQ.0) GO  TO  20 

READ  (5, FORI)  (DATA  (L)  ,L  =  1,II) 

20  IM=I 1+  1 

READ  (5 , FOR2 )  (DATA  (L)  ,L  =  IM,NCAS) 

C 

C  FIND  THE  VARIANCE,  MEAN  AND  STANDARD  DEVIATION  OF  THE  DATA. 

C 

CALL  VARSD  (NCAS, DATA, EM, S,SD) 

C 

C  PRINT  OUT  THE  JOB  PARTICULARS. 

C 

KQ=I ABS (KP) 

WRITE  (6, 103) JOB, NCAS, KQ, SIGN 
WRITE(6, 104) EM , S , SD 
C 

C  PRINT  OUT  DATA  IF  REQUIRED. 

C 

IF  (KP. EQ. 0) GO  TO  40 
MN=  1 

25  NN=MN+ 1  1 

IF (NN.GT. KQ) NN=KQ 
WRITE(6,105)  (J,J=MN,NN) 

WRITE  (6 , 106)  (DATA  (I)  ,I  =  MN, NN) 

MN=M  N  + 1 2 

IF (MN.LE. KQ) GO  TO  25 
IF (KP.LT. 0) GO  TO  40 
C 

C  PLOT  THE  DATA. 

C  THE  FIRST  (LB)  DATA  VALUES  ARE  SET  TO  VALUES  WHICH  ARE  HOPEFULLY 

C  BEYOND  THE  RANGE  OF  THE  PLOT. 

C 

IF  (LB. LT. 0) LB=0 
DO  30  1=1, NCAS 
30  PLOT (I+LB) =DATA (I) 

IF  (LB.LE.0) GO  TO  38 
DO  35  1=1 , LB 
35  PLOT  (I) =1000000.0 
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38  NTOT=NCAS+LB 

CALL  PLTDAT  (PLOT, NTOT) 

C 

C  RFMOVE  THE  MEAN  FROM  THE  DATA  FOR  OSE  BY  THE  F.F.T.  AND  STORE 

C  IN  ARRAY  V, 

C  NORMALIZE  THE  DATA  ARRAY  FOP  THE  M.E.M.  AND  STORE  IN  ARRAY  W. 

C 

UO  DO  45  1=1, NCAS 
V  (I)  =DAT A  (I)  -EM 
u5  W  (I)  =V  (I)  /SD 
C 

C  NOW  PERFORM  THE  ANALYSIS  REQUIRED  AND  THE  PLOTTING. 

C  IF  BOTH  FFT  AND  MEM  OUTPUTS  ARE  PLOTTED,  SPECIFY  PARAMETERS  FOR 

C  THE  FFT  PLOT  FIRST.  IF  BOTH  UNSMOOTHED  AND  SMOOTHED  POWER 

C  SPECTRUM  FROM  THE  FFT  ARE  PLOTTED,  SPECIFY  PARAMETERS  FOP  THE 

C  UNSMOOTHED  SPECTRUM  FIRST. 

C 

KT2=KT+u 
NZ  P=  NZERO/2+ 1 

GO  TO  (50, 55, 60,65,  70, 75,80)  , KT2 
50  CALL  FASTFT (NCAS , NX, NZEPC, NOM, NFIL, SIGN, DT, V) 

CALL  BURG ME (NCAS, DT , W, MP, MEMN, ME MX) 

GO  TO  90 

55  CALL  BURGME  (NCAS , DT, W, MP, MEMN, MEMX) 

GO  TO  90 

60  CALL  FASTFT  (NCAS , NX, NZERO, NOM, NFIL, SIGN, DT, V) 

GO  TO  90 
65  WPITE (6 , 107) 

GO  TO  90 

70  CALL  FASTFT (NCAS , NX, NZERC, NOM, NFIL, SIGN, DT, V) 

CALL  PLTFFT (NZP , KN) 

GO  TO  90 

75  CALI,  BURGME  (NCAS,  DT,  W,  MP,  MEMN,  MEMX) 

CALL  PLTMEM ( MP) 

GO  TO  90 

80  CALL  FASTFT  (NCAS , NX, NZERO, NOM, NFIL, SIGN, DT, V) 

CALL  PLTFFT (NZP, KN) 

CALL  BURGME  (NCAS, DT, W, MP, MEMN,  MEMX) 

CALL  PLTMEM (MP) 

90  CONTINUE 
GO  TO  5 
95  STOP 
END 


on  non 


Table  B.2.  Subroutine  VARSD,  as  described  in  Chapter  6 


SUBROUTINE  VARSD ( NCA S , DAT  A, E M, S r S D) 

*****#$*'*******♦*****#***$*#***  •fr***#*************#**:**#***^***:**** 

FIND  THE  VARIANCE  AND  STANDARD  DEVIATION. 

************  ************************************************  ****** 

110  FORMAT(1H1,20X, ’ *****  ZERO  OR  NEGATIVE  VARIANCE  *****’) 

*  ***  ***  *  ****  ****  *  **********************************  ******  ********* 

DIMENSION  DATA  (NCAS) 

RC  A=  FLOAT  ( NCAS) 

SSQ=0 . 0 
EMN=0 . 0 

DO  10  1= 1 , NC AS 
EMN=EMN+D AT A (I) 

10  SSQ=SSQ  +  DATA  (I)  *DATA  (I) 

EM=EMN/RC  A 
S=SSQ/RCA-EM*EM 
IF  (S.GT.O.) GO  TO  20 
S=0. 0 

WRITE  (6,110) 

SD=0 . 0 
RETURN 

20  SD  =  SQRT  (3) 

RETURN 

END 
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Table  B.3.  Subroutine  FASTFT,  as  described  in  Chapter  4 . 


c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  FASTFT ( NC AS , NX , NZERO , N CM, NFIL, SIGN  f  DT , DATA) 

X  IS  AN  ARRAY  OF  COMPLEX  NUMBERS. 

A  LABELLED  COMMON  BLOCK  WILL  HOLD  THE  RESULTS  FOR  PLOTTING. 


ADD  ZEROS  AT  THE  BEGINNING  AND/OP  THE  END  OF  THE  DATA  TO  OBTAIN 
THE  REQUIRED  LENGTH  NZERO. 

FIND  THE  NEW  MEAN  AND  STANDARD  DEVIATION,  AND  NORMALIZE  THE  DATA. 

119  F0RMAT(1H1, '  *****  NZERO  TOO  SMALL  ***♦*') 

120  FORMAT ( 1  HI FOP  THE  F.E.T.,  THE  MEAN  =  ',E10.6,'  TOTAL  VARIANCE  = 

* ' r  FI  3 . 5 , '  STANDARD  DEVIATION  =  • , F 1 0 . 4 r / ,  *  AFTER  APPLICATION  OF  A 
* • , iu , «  POINT  TAPER. » ) 

121  FORMAT ( 1H-, 12X, ' TIME’ , 8X,  ' REAL  X*  , 5X , * IMAGINA RY  X*,/) 

122  FOFM  AT (  1  H  , 1 IX, 14,2 ( 7X , F8.3)  ) 

123  FO  RM  AT  (  1 H- , '  ALL  SMOOTHING  APPLIED  USING  A', 13,'  POINT  DANIELL  (BO 

♦XCAR)  FILTER ’,//, 58X  ,» UNSMOOTHED  SMOOTHED  SMOOTHED  W 

WEIGHTED  SMOOTHED *,/, •  FREQUENCY  COS INE •  ,  1 0 X ,  1 S I NE • , 1 0 X ,  • 

♦PHASE', 8X, 'POWER' ,10Xr ' PO WER • , 8 X , ’ A MPL .  SQ.  FREQUENCY  PHASE 

*',//) 

12  4  FORM  AT  ( 1 H  , F9 . 5 , 1 X , 2 El  5. 5 , F 1 2 . 3 , 3 E 1 5 . 5 , F 1  1 . 5 , F 1 2 . 3) 

125  FORMAT  ( 1  HO  , *  THE  SUMS  ARE ' , 39X, 3E 1 5 . 5, //, '  ALSO  NOTE  THAT  THE  FIR 
*ST  AND  LAST' ,13,’  POINT(S)  OF  THE  SMOOTHED  POWER,  AMPL.  SQ.,  PHASE 
*',/,  ’  AND  WEIGHTED  FREQUENCY  COLUMNS  HAVE  NOT  BEEN  SMOOTHED’) 

126  EO  RM  AT  ( 1 H-  ,  '  ALL  SMOOTHING  APPLIED  USING  A’  ,13,'  POINT  DANIELL  (BO 

*XC AR)  FILTER ',//, 58X ,» UNSMOOTHED  UNSMOOTHED  WEIGHTED  SMO 

♦OTHED ' , '  FREQUENCY  COSINE ’, 1 0 X ,’ S INE ’, 1 0 X ,' PHA S E 8X ,’ POW E 

*R' ,8X, ' AMPL.  SQ.  FREQUENCY  PHASE',//) 

127  FORMAT  ( 1 H  , F9 . 5 , 1 X , 2E 1 5. 5 , F 1 2. 3 , 2E 1 5 . 5 , F 1 1 . 5, F 1 2 . 3) 

128  FORM  AT  ( 1  HO  ,  ’  THE  SUMS  APE ' , 3 9X , 2E 1 5 . 5 , //,  '  ALSO  NOTE  THAT  THE  FIR 
*ST  AND  LAST', 13,’  POINT(S)  OF  THE  SMOOTHED  PHASE',/,'  AND  WEIGHTE 
* D  FREQUENCY  COLUMNS  HAVE  NOT  BEEN  SMOOTHED') 

129  FORM  AT  (1 H-r '  NO  SMOOTHING  HAS  BEEN  A PPLI E D ' , // , 1 0 X , ' FREQUE NC Y ' , 7 X  , 
♦  •COSINE' r12X,'SINE* ,12X, ' PHASE' , 10X, ' POWER' , 1  OX, ' AMPL.  SQ.  ' ,//) 

130  FORMAT  ( 1 H  , 8 X , E9 . 5 , 2 E 1 7. 5 , F 1  a. 3 , 2E 1 7 . 5 ) 

131  FORM  AT  ( 1  HO , ’  THE  SUMS  A  RE  * , 5 2X , 2E 1 7 . 5) 


DIMENSION  CX  (2  0  5  0)  , SX  (20  50)  , PHZ  (2Q5  0)  , DATA (NCAS)  ,XX(U096)  ,XC(2050) 
*, XS(2050)  , WOP (20  5  0)  ,ZHP(2  050) ,QERF(20  50)  , AM  PS Q  (2  0  50) 

COMPLEX  X  (4  096) 

COMMON  /OUTPLT/  FREQ  (2 050) , POW  (2 0 50 ), SO W  (2050 ) 

TAPER  THE  DATA  AT  BOTH  ENDS  (IF  DESIRED)  USING  A  COSINE  BELL. 

IF (NOM. LE. 0) GO  TO  16 
CALL  TAPER (NCAS, DATA , NCM) 

CALCULATE  THE  MEAN  AND  VARIANCE  AND  ASSIGN  THE  DATA  TO  THE  COMPLEX 
C  ARRAY  X. 

C 

16  DO  17  1=1, NCAS 

17  XX  (I) =DATA  (I) 

IF (NZERO. EQ. NCAS) GO  TO  30 
IF  (NZERO. GT.  NCAS)  GO  TO  20 


•  .w. 
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WRITE (6,119) 

RETURN 

20  MNM=NCAS+1 

DO  25  I=MNM, NZERO 
25  XX(I)=0.0 

30  CALL  VARSD  (NZERO, XX, EMM, SS,  SD) 

WRITE (6 ,120} EMM, SS,SD, NOM 
DO  35  1=1, NZERO 
35  X  (I)  =  (XX  (I) -EMM) /SD 

CALL  NLOGN  (NX, NZERO, X, SION) 

WE  NOW  HAVE  THE  FOURIER  TRANSFORM  OF  THE  DATA  STORED  IN  ARRAY  X. 

IF  (SIGN. LT. 0.0) GO  TO  50 
WRITE  (6 ,121) 

DO  45  1=1, NZERO 
45  WRITE  (6, 122)  I, X (I) 

RETURN 

50  RNZ  =  FLOAT (NZERO)  *DT 
NZP=NZERO/2+1 
DO  60  1=1, NZP 
CX  (I)  =REAL  (X  (I)  ) 

SX  (I)  =  A I M  AG  (X  (I)  ) 

60  FREQ  (I) =FLOAT (1-1) /RNZ 

CALCULATE  THE  POWER  AND  PHASE. 

CALL  POLAR  (NZP, CX,SX ,P0W, PHZ) 

FAM=2.0*SS*FLOAT  (  NZERO) /FLOAT (N CAS) 

APPLY  DANIEL!  FILTER  WHERE  REQUIRED. 

A  MINIMUM  VALUE  OF  .000000011  IS  ASSIGNED  TO  THE  POWER  (AFTER 
PRINTING)  IN  ORDER  TO  PREVENT  ERRORS  IN  THE  PLOTTING 
PROCEDURE. 


BX=. 11E-08 
PS  UM  =0 . 0 
SSUM=0. 0 
ASUM=0 . 0 

IF  (NFIL. EQ. 0) GO  TO  80 
MFIL=I ABS (NFIL) 

KA=  (MFIL-1) /2 

CALL  DAN  I  (MFIL ,NZP,CX, XC) 

CALL  DANI  (MFIL , NZP, SX, XS) 

CALL  POLAR  (NZP, XC , XS , WCP, ZH P) 

CALL  DANY  (MFIL, NZP, FREQ, POW,QERF) 

IF  (NFIL.GT. 0) GO  TO  70 
CALL  DANI  (MFIL , NZP, POW, SOW) 

THE  RESULTS  ARE  PRINTED. 

WRITE  (6 , 123)  MFIL 
DO  66  1=1, NZP 
AMPSQ(I) =FAM*SOW (I) 

WPITE(6,12U)  FREQ  (I)  ,CX(I)  ,SX(I)  , PH  Z  (I )  ,POW(I)  , SOW (I)  , AMPSQ (I)  ,QERF 
*(I),ZHP(I) 

PSUM=PSUM+POW (I) 

SSUM=SSUM+SO W (I) 

ASUM  =  AS  UM  +  AM PSQ (I) 

IF  (POW (I)  . LT.BX)  POW  (I)  =3X 
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IF  (SOW (I)  .LT.BX)  SOW  (I) =BX 
CONTINUE 

WRITE  (6 ,125) PSUM,SSUK,ASUM,KA 
RETORN 

WRITE  (6, 126)  MFIL 
DO  75  1=1, NZP 
AMPSQ  (1) =FAM*POW (I) 

WRITE (6, 12?)  FREQ (I)  ,CX  (I)  ,SX  (I) 

♦  (I) 

PSUM=PSUM+POW (I) 

ASUM=ASUM+AMPSQ (I) 

IF  (POW (I)  .LT.BX)  POW (I) =BX 
CONTINUE 

WRITE  (6, 12  8)  PSUM, ASUM, KA 
RETURN 

WRITE (6, 1 29) 

DO  85  1  =  1 , NZ P 
AMPSQ (I) =FAM*POW (I) 

WRITE  (6,  130)  FREQ  (I)  ,CX(I)  ,SX  (I) 
PSUM=PSUM+POW {I) 

ASUM=ASUM+AMPSQ (I) 

IF  (POW(I)  .LT.BX)  POW 1 1)  =BX 
CONTINUE 

WRITE  (6, 131)  PSUM, ASUM 

RETURN 

END 


,  PHZ  (I)  ,POW(I)  ,  AMPSQ  (I)  ,QERF(I)  ,ZHP 


,  PH  Z  (I)  ,POW(I)  ,  AMPSQ(I) 
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Table  B.4.  Subroutine  DANI,  as  described  in  Chapter  u. 


SUBROUTINE  DAN I (NFIL, N, POW, BMP) 

C  ******^*****************  ******************************  **#*#«:***** 

C  NFIL  =  NUMBER  OF  FILTER  POINTS:  3,  5,  7,  9  OR  11. 

C  N  =  NZEFO/2  +  1. 

C  POW  =  UNFILTERED  ARRAY  INPUT. 

C  BMP  =  FILTERED  ARRAY  OUTPUT. 

C 

C  NOTE  THAT  THE  FIRST  AND  LAST  (NFIL-1)/2  POINTS  OF  THE  ARRAY  WILL 

C  EMERGE  UNFILTERED. 

Q  **  #  **  *****  ***  *****  *  *******  ****  *  ***********************  ********** 

400  FORMAT (1H1 ,» *****  NUMBER  OF  FILTER  POINTS  EXCEEDS  11  **♦**') 

****************************************************************** 

DIMENSION  POW  (N)  , BMP  (N) 

IF  (NFIL. GT. 1 1) GO  TO  50 
DO  10  M= 1 r  N 
10  BMP(M) =0.0 

FII.=  FLOAT  (NFIL) 

KA= ( NFIL- 1 ) /2 
NA=K A+ 1 
NB=N -K  A 
DO  30  I  =  N  A, NB 
NC=I -K A 
ND=I +K A 
DO  20  J  =  NC , N  D 
20  BMP(I) =BMP (I) +POW (J) 

30  BMF(I) =EMP (I) /FIL 
NK=N+1 

DO  40  1=1, KA 
BMP(I)  =POW  (I) 

40  BMP(NK-I) =POW (NK-I) 

RETURN 

50  WRITE  (6,4  00) 

RETURN 
END 


* 
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Table  B.5.  Subroutine  DANY,  as  described  in  Chapter  u. 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  DANY (NFIL, N,FR,PO,QE) 

MODIFIED  CANIELL  FILTER  FOR  RESOLVING  THE  PEAK  FREQUENCIES  USING 
THE  POWER  AS  A  WEIGHTING  FACTOR. 

NFIL  =  NUMBER  OF  FILTER  POINTS:  3,  5,  7,  9  OR  11. 

N  =  NZERO/2  +  1. 

FR  =  ARRAY  TO  BE  FILTERED. 

PO  =  ARRAY  USED  AS  WEIGHTS. 

QE  =  FILTERED  ARRAY  OUTPUT. 

NOTE  THAT  THE  FIRST  AND  LAST  (NFIL-1)/2  POINTS  OF  THE  ARRAY  WILL 
EMERGE  UNFILTERED. 

450  FORMAT  (1H1  , •  *****  NUMBER  CF  FILTER  POINTS  EXCEEDS  11  ******) 

** ***************************** **********  ***************  ********** 


DIMENSION  FR  (N)  ,  PO  (N)  ,  QE  (N) 
IF  (NFIL. GT. 11) GO  TO  50 
DO  10  M  = 1 ,  N 
10  QE (M ) =0 . 0 

FIL= FLOAT (NFIL) 

KA  = ( N  FIL - 1 ) /2 

NA=KA+ 1 

NB=N -K A 

DO  30  I  =  NA,  NB 

NC=I -K A 

ND=I+KA 

POW=0. 0 

DO  20  J=NC,ND 

QE  (I)  =QE  (I)  +FR  (J)  *PO  (J) 

20  POW=  POW+  PO ( J) 

30  QE  (I)  =QE  (I) /POW 
NK=N+ 1 

DO  40  1=1 ,  K  A 
QE  (I)  =FR  (I) 

40  QE  (NK-I) =FP (NK-I) 

RETURN 

50  WRITE (6 , 450) 

RETURN 

END 
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Table  B.6.  Subroutine  TAPER,  as  described  in  Chapter  4. 


SUBROUTINE  T  APER  ( N ,  X  ,  LC) 

C  S************************************** #$**** **$******#*********** 

C  N  =  TOTAL  SIZE  OF  ARRAY  X. 

C  X  =  DATA  ARRAY. 

C  LC  =  NUMBER  OF  FILTER  POINTS  IN  H  AL  p  THE  FILTER,  MAXIMUM  100. 

C 

DIMENSION  X  (N)  , T  1 1 00) 

PI  =  3. 141592654 
CH  =  FLOAT  (LC- 1 ) 

DO  10  1=1, LC 

10  T  (I)  =  (COS  (PI *FLOAT  (I -LC)  /CM)  +1 .  0)  /2. 0 

T ( 1 )  =  0.0  ,  T  ( LC )  =  1.0 

DO  20  1=1, LC 
K=N-I+ 1 

X  (I)  =X  (I)  *T  (I) 

20  X  (K)  =X  (K)  *T  (I) 

PETURN 
END 
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Table  B.7.  Subroutine  POLAR,  as  described  in  Chapter  4. 


SOB ROUTINE  POLAR (L, RE, XIM, AMP, PHZ) 

*******$:*******#**^**:(c*:fr**#*:«c:*****$#:***:Cc**;<::(c*:ft***:C£**:«r************ 

OBTAINS  THE  PH  AS  F  AND  THE  POWER  (AMPLITUDE  SQUARED)  . 
#*$#**#*****$******#*#*****#;*#*##$*$********#**##*:*#;****♦:*«#*  +  **** 

DIMENSION  RE  (L)  ,  XIM  (L)  , AM P (L)  , PHZ  (L) 

PI=3 . 141592654 
ANG= 180. /PI 
DO  120  1=1, L 

AMP  (I) =2,0*  (RE (I) *RE (I) +  XIM (I)  *XIM (I)  ) 

IF  (XIM (I ) )  10,20,30 

10  IF  (RE  (I) )  40,50,60 
20  IF  (PE (I) )  70,80,60 

30  IF  (PE  (I) )  100,90,60 

4  0  PHZ  (I)  =  AT  AN  (XIM  (I ) /P  E  ( I)  )  -  PI 
GO  TO  110 
50  PHZ (I) =-90.0 
GO  TO  120 

60  PHZ (I) =ATAN  (XIM (I)/PE  (I) ) 

GO  TO  110 
70  PHZ(I) =-180. 0 
GO  TO  120 
80  PHZ(I) =0.0 
GO  TO  120 
90  PHZ (I) =90.0 
GO  TO  120 

100  PHZ  (I) = AT AN  (XIM (I) /RE (I)  )  +  PI 
110  PHZ (I) =PHZ (I) *ANG 
120  CONTINUE 
RETURN 
END 
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Table  B.8.  Subroutine  NLOGN,  as  described  in  Chapter  4 


SUBROUTINE  NLOGN ( N, LX, X, SIGN) 

C  ****************************************************************** 

C  THE  RETORN  OOTPUT  IN  ARRAY  X  IS  THE  FOURIER  TRANSFORM  OF  THE  DATA. 

C  M  IS  DIMENSIONED  TO  HAVE  AS  MANY  VALUES  AS  THE  LARGEST  VALUE  OF  N. 

C  (FROM  ROBINSON,  1967) 

Q  ************************************************* ******* ********** 

c 

DIMENSION  M  1 1 2) 

COMPLEX  X  (LX)  , WK , HOLD, G,C MPLX 
DO  10  1=1  ,N 
10  M  (I)  =2**  (N-I) 

FLX=FLOAT  (LX) 

PIX=SIGN*2. *3. 141 59265U/FLX 
DO  40  L=  1  ,  N 
NBL=2**  (L-1) 

LBL=LX/N  BL 

LBH=LBL/2 

K=0 

DO  40  I BL= 1 , NBL 
FK=K 

V= PI X*FK 

WK=C  M  ILX  (COS  (V)  ,  SIN  (V)  ) 

IST=LBL*  (IBL-1) 

DO  20  1=1, LBH 
J=IST+I 
JH=J+LBH 
Q=X  (JH) *WK 
X (JH) =X  (J) -Q 
20  X  (J)  =X  (J)  +Q 
DO  30  1=2, N 
II=I 

IF  (K.LT. M  (I)  ) GO  TO  40 
30  K=K-M(I) 

40  K=K+  M  (II) 

K=  0 

DO  70  J= 1 , LX 
IF  (K.LT. J) GO  TO  50 
HOLD  =  X  (J) 

X  ( J)  =X  (K+  1 ) 

X  (K+ 1) =HOLD 
50  DO  60  1=1, N 
II  =  I 

IF  (K.LT. M  (I)  ) GO  TO  70 
60  K=K- M (I ) 

70  K=K  +  M  (II) 

IF(SIGN.GT.O.O) RETURN 
DO  80  1=1, LX 
80  X  (I)  =X  (I)  /FLX 
RETURN 
END 
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Table  B.9.  Subroutine  BURGME,  as  described  in  Chapter  5 


SUBROUTINE  BURGME  (N, DT , DATA , NZ P , MEMN , MEMX) 

C  USING  THE  MAXIMUM  ENTROPY  METHOD,  WE  HOW  CALCULATE  A  DIFFERENT  SET 

C  OF  SPECTRAL  ESTIMATES. 

C  THE  AKAIKE  FINAL  PREDICTION  ERROR  IS  USED  TO  DETERMINE  THE  OPTIMUM 

C  NUMBER  OF  FILTER  COEFFICIENTS.  I?  THE  NUMBER  OF  COEFFICIENTS 

C  IS  SPECIFIED,  THE  AKAIKE  FPE  IS  BYPASSED,  AND  FPE l M)  AND 

C  RFPE(M)  ARE  SET  AS  ZERO. 

C 

C  N  =  NUMBER  OF  DATA  POINTS. 

C  DT  =  TIME  INTERVAL  CHOSEN  AS  BEFORE. 

C  DATA  =  NORMALIZED  DATA. 

C  NZ P  =  NUMBER  OF  FREQUENCY  AND  POWER  VALUES  CALCULATED. 

C  MEMN  =  NUMBER  OF  FILTER  COEFFICIENTS  IF  NOT  CHOSEN  AUTOMATICALLY 

C  BY  AKAIKE’ S  FPE,  MAXIMUM  (N-2). 

C  ME MX  =  EXPANSION  FACTOR  FOR  LOWER  END  OF  FREQUENCY  SPECTRUM. 

C 

C  A  LABELLED  COMMON  BLOCK  WILL  HOLD  THE  RESULTS  FOR  PLOTTING. 

Q  ^cjc*#*:******#**************##*#*:#********  ****$*****$******:****♦**** 

145  FORMAT  (1H1, ’THE  NUMBER  OF  FILTER  COEFFICIENTS  CALCULATED  IS’,I5,’. 

*  THE  EXPANSION  FACTOR  IS’,15,’.  DT  =  ’,F5.3) 

146  FORMAT (1H1 , ’THE  NUMBER  OF  FILTER  COEFFICIENTS  WILL  BE  CHOSEN  AUTOM 
* ATIC ALL Y .  THE  EXPANSION  FACTOR  IS’,15,’.  DT  =  *,F5.3) 

147  FORM  AT (1H-, ’ THE  PROCESS  OF  CONVERGENCE  TO  MINIMUM  FPE  NOW  FOLLOWS’ 

*) 

148  FORM  AT (1H0, 4X, ’M  FPE  ( M)  P(M)*,/) 

149  FORMAT  1 1 H  ,I5,2E15.5) 

150  FORMAT (1H1 ,’ FOR  M  =  ’,15,’  WE  HAVE  FPE (M)  =  *,E15.5,’  AND  P*M)  =  • 
*, E15. 5,/, ’THE  FILTER  COEFFICIENTS  ARE  ’) 

151  FOPMATlIH  ,  6F1  0 . 5) 

152  FORMAT  (1H1 ,’ FPE  (0)  =  ’,E15.5,»  FPE  ( M)  =  *,E15.5,’  RELATIVE  VALUE  ( 
*RFPE  iM)  )  =  ’ f E15. 5,/) 

1 53  FORMAT (1H0, 12 X, ’ FREQUENCY’ , 1 IX, ’ POWER’ ,/) 

154  FO  RM  AT ( 1 H  , F 20 . 5 , E2 0 . 5 ) 

****$#*********#*****#*********##**#*$ ♦^S***#*###*'***#*****#****** 

DIMENSION  DATA  ( N )  ,  Y  A  (4096),YB(4096),AA(4096)  ,  AB  (  409  6)  ,CA  <  4  096) 
COMMON  /MEMOUT/  FR EQ  (4  G97 )  ,  P  (U  097  ) 

INITIALIZE  THE  PROCESS. 

MT=2  5 

RCA=FLOAT (N ) 

IS  =  0 
SUM=0 . 0 
DO  15  1=1, N 

15  SUM=SUM  +  DATA (I)  *DATA  (I) 

PM=SUM/RCA 

FFEO=PM*FLO AT (N + 1 ) /FLO AT l N- 1 ) 

NMl=N-1 
YA  1 1 ) =DATA (1 ) 

YB  (NM  1 )  =DATA  (1) 

DO  2  0  1=  2 , N  M 1 
YA  (I )  =  D ATA  (I) 

20  YB  ll -  1 )  =DATA  (I) 


* 


^2G 


c 

C  WE  WISH  TO  FIND  THE  MINIMUM  VALUE  OF  THE  AKAI KE  FPE  AND  THE 

c  COEFFICIENTS  ASSOCIATED  WITH  THIS  MINIMUM,  UNLESS  THE  NUMBER 

C  OF  FILTER  COEFFICIENTS  HAS  BEFN  SPECIFIED. 

C 

IF  lMEMN.EQ.O)GO  TO  25 
WRITE  (6 , 145)  MEMN, MEM X , DT 
GO  TO  30 

25  WRITE  (6, 1 46)  MEMX, DT 
WRITER,  147) 

WRITE  (6,148) 

30  DO  75  M= 1 , N 
M A=M - 1 
MB  =  N-M 

IF  (M .EQ. 1) GO  TO  45 
DO  35  1=1, MA 
35  AB  (I)  =A A  (I) 

DO  40  1=1, MB 

YA  (I)  =  YA(I)  -  AB  (M-1)  *YB  (I) 

40  YB (I) =YB  (I+1J-AB  (M-1) *  YA  (1+1) 

45  DOM  =  0 . 0 
DEN=0 , 0 
DO  50  1=1, MB 
DOM=DOM  +  Y A  (I ) *YB (I) 

50  DEN  =  DEN+  YA  (I)  *  YA  (I)  +YB  (I)  *YB  (I) 

AA  (M)  =  2 . 0  *  DO  M/ D E N 
PM=PM*(1.0-AA(M)  *  A  A  (M)  ) 

IF  (M .  EQ. 1) GO  TO  70 
DO  55  1=1, MA 

55  AA  (I)=AB(I)-AA(M)  *A3  (M-I) 

IF  (MEMN-M)  57,61,75 
C 

C  THE  RESIDUAL  SUM  OF  SQUARES  IS  NOW  CALCULATED  IF  THE  AKAIKE  FPE  IS 

C  TO  BE  DETERMINED  BY  THE  SUBROUTINE. 

C 

57  CALL  SMSQR (D  ATA,  N  ,  AA  ,  M  ,  SM) 

XFPE=SM*?LOAT (N+M+ 1) /FLOAT ( N-M- 1 ) 

IF  1XFPE.LT. FPEM) GO  TO  60 
IS=IS+1 
GO  TO  73 
C 

C  STORE  THE  LATEST  SET  OF  VALUES. 

C 

60  FPEM=XFPE 

61  PN=PM 
MM  =  M 

DO  62  1=1, M 

62  CA  (I) =AA  (I) 

IF  (M.EQ. MEMN) GO  TO  80 
IS  =  0 

GO  TO  73 
C 

C  CALCULATE  THE  INITIAL  VALUES. 

C 

70  IF  (MEMN.GT.O)GO  TO  71 
CALL  SMSQR(DATA,N,AA,M,SM) 

FPEM=SM*FLOAT (N+ 2) /FLOAT (N-2) 

71  MM= 1 
PN  =  PM 

CA  ( 1 )  =  A A  {"l) 


*  *  )  , 


IF  (MEMN.EQ. M) GO  TO  3  0 
IF  (MEMN.GT.O)GO  TO  75 
XFPE=FPEM 

73  WRITE  (6, 149)  M, XFPE, PM 
IF  (I S . GT. MT) GO  TO  80 
75  CONTINUE 
C 

C  NOW  WE  CALCULATE  THE  FREQUENCY  AND  POWER  VALUES. 

C  A  MINIMUM  VALUE  OF  .000011  IS  ASSIGNED  TO  THE  POWER  (AFTEP 

C  PRINTING)  IN  ORDER  TO  SIMPLIFY  THE  CHOICE  OF  PLOTTING 

C  PARAMETERS. 

C 

80  IF  (MEMN.GT. 0) FPEM  =  0. 0 
WRITE  (6  , 150)  MM , FPEM , PN 
WRITE (6,151)  (CA  (I)  ,1  =  1,  MM) 

T2PI=DT*6. 2831 85307 
BX=. 000011 

RNZ=2.0*DT*FLOAT (N) * FLOAT (ME MX) 

PT=PN*DT 
NZP=N+ 1 

R  FPE=FPEM/FPEO 

WRITE  (6,152)  FPEO, FPEM, RFPE 

WRITE (6 , 153) 

DO  90  1=1 ,NZP 

FREQ  (I) =FLOAT (1-1 ) /RNZ 

ARG=T2PI*FREQ (I) 

SREAL= 1 .0 
SIMAG=0 .0 
DO  85  J  = 1 , M M 

SREAL=5RE AL-CA ( J) *COS(ARG*FLOAT ( J)  ) 

85  SIMAG  =  SIM AG  +  CA (J) *SIN (ARG*FLOAT (J)  ) 

P  (I) =PT/  (SREAL*SREAL+SIMAG*SIMAG) 

WRITE  (6, 15  4)  FREQ (I)  , P (I) 

IF  (P  (I)  .LT. BX) P (I) =BX 
90  CONTINUE 
RETURN 
END 
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Table  B. 10.  Subroutine  SMSQR,  as  described  in  Chapter  5. 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  S MS QR { X , N , A , M, SM) 

THIS  PROGRAM  CALCULATES  A  RESIDUAL  SUM  OF  SQUARES  NEEDED  TO 
COMPUTE  THE  AKAIKE  FEE. 

X  =  DATA  ARRAY. 

N  =  NUMBER  OF  DATA  VALUES. 

A  =  COEFFICIENT  ARRAY. 

fi  =  NUMBER  OF  COEFFICIENT  VALUES. 

SM  =  RESIDUAL  SUM  OF  SQUARES. 

♦  *$***  +  ****:***:$:*  ************  ***■*  *******$******♦************#* 


10 

20 

30 


DIMENSION  X  (N)  , A (M) 

SUM=0 . 0 

DO  30  1=1, N 

SUN=0 . 0 

DO  10  J=1,M 

IF  (J .GE. I) GO  TO  20 

SUN=  SUN+  A  {J)  *X  A I  —  *1 ) 

CONTINUE 

SUM=SUM+(XlI)-SUN)**2 
SM= SUM/FLOAT (N) 

RETURN 

END 


p-  4 


o  n  n  no 
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Table  B«11.  Subroutine  PLTDAT,  as  described  in  Chapter  6. 


SUBROUTINE  P LTDAT  (DATA , N) 

C  *******<'*««*^*^****«****^*********^* ♦*♦♦***«****♦*«♦**** ********** 

C  THE  RAW  DATA  IS  STORED  IN  ARRAY  DATA. 

C  ARRAY  YRS  IS  CALCULATED  AND  CONTAINS  VALUES  FOR  THE  ABSCISSA  AXIS. 

C  ALPHA  CONTAINS  THE  TITLE  AND  AXIS  LABELS. 

C  Z  PLOTTING  FACILITY  NOT  AVAILABLE. 

C  PARAMETERS  INPUT  IN  THIS  SUBPROGRAM  (FROM  THE  FILE  ASSIGNED  TO 

C  UNIT  7)  : 

C  YR1  =  FIRST  YEAR  TO  BE  PLOTTED  (WITH  DECIMAL  PERIOD  IF  FIRST 

C  MONTH  IS  NOT  JANUARY). 

C  YRIN  =  INCREMENT  OF  TIME  FOR  EACH  DATA  POINT. 

C  NF  =  AS  IN  CSLIB  (COMPUTING  SCIENCE  LIBRARY  WRITE-UP  FOR  CGPL) . 

C  KA  AND  KB  AS  IN  CSLIB  EXCEPT  FOR  E A=KB= 1 ,  WHERE  -1  ROTATES  THE 

C  NUMBERS  90  DEGREES  W.R.T.  THE  AXIS,  FOR  ARITHMETIC  AXIS  ONLY. 

C  KC  AS  IN  CSLIB  FOR  KC  =  3  TO  6  (ALSO  SEE  COMMENTS  IN  OTHER 

C  PLOTTING  SUBROUTINES). 

C  FOR  KC  =  1  OR  2  SYMBOL  HEIGHT  DETERMINED  BY  HT, 

C  FOR  KC  =  7  SYMBOLS  OF  HEIGHT  RT  APE  PLOTTED  WITH  NO  LINE, 

C  ACCORDING  TO  SYD  AND  SYS. 

C  FOR  KC  =  8  SYMBOLS  OF  HEIGHT  HT  ARE  PLOTTED  AND  JOINED  BY 

C  STRAIGHT  LINE  SEGMENTS  ACCORDING  TO  SYD  AND  SYS. 

C  HA, HB , HC  AS  IN  CSLIB. 

C  VA,VB,VC  AS  IN  CSLIB. 

C  HT  HEIGHT  OF  PLOTTED  SYMBOLS  IN  DECIMAL  INCHES 

C  (DUMMY  VARIABLE  USED  WHEN  KC  =  3  TO  6) . 

C  SYD  DELAY  (IN  DECIMAL  INCHES)  BEFORE  FIPST  SYMBOL  IS  PLOTTED 

C  ON  THE  LINE. 

C  SYS  SEPARATION  (IN  DECIMAL  INCHES)  BETWEEN  SYMBOLS  PLOTTED 

C  ON  THE  LINE. 

C  NOTE:  SYD  AND  SYS  MAY  BE  DUMMY  VAPIABLES  WHEN  KC  =  1  TO  6. 

C  ****************************************************************** 

600  FORMAT  (F9.4,E6.4) 

601  FO  RM  AT ( 2  0  A4) 

602  FORMAT  (UI5,6F10. 3) 

603  FORMAT (3F10. 3) 

****************************************************************** 

DIMENSION  DATA  (N)  , YRS  (4096)  , ALPHA (2  0) 

READ  (7,  600)  YP1 , YRIN 
DO  10  1=1, N 

10  YRS  (I) =YR 1+  (I- 1) *YRIN 

READ  (7, 601)  (ALPHA  (I)  ,1  =  1,20) 

READ  (7,6  02) NF, KA , KB , KC , HA , HB , H C , V A , VB , VC 
READ  (7,603) HT, SYD, SYS 

CALL  THE  PLOTTING  SUBROUTINE. 

CALL  PLTG  (YRS,DATA,DATA,N,NF,KA,KB,KC, KC , HA , HB , HC , V A , VB , VC , ALPH A , 6 
*,HT,SYD,SYS, 1) 

RETURN 
END 


ono  non  on 


Table  B.12.  Subroutine  PLTFFT,  as  described  in  Chapter  6. 


SUBROUTINE  PLTFFT (N ,  M) 

c  a***************^**********  *******  *********•.*#*:#***:****#:((  *********** 

C  THE  PARAMETERS  TO  BE  READ  IN  THIS  SUBROUTINE  (FRGM  UNIT  7)  ARE  THE 

C  SAME  AS  IN  PLTDAT. 

C  IF  NO  PLOT  IS  WANTED  AT  THIS  TIME  BUT  THE  INFORMATION  IS  TO  BE 

C  STORED,  SET  KC  NEGATIVE  AND  ASSIGN  A  FILE  TO  UNIT  1. 

Q  *******-/********************************************************** 

6  10  FORMAT  (20 A4) 

611  FORMAT (4l5,6?10.3) 

612  FORMAT (3F10. 3) 

613  FORMAT (8F10. 5) 

*******************  *********************************************** 
DIMENSION  ALPHA  (20) 

COMMON  /OUTPLT/  FREQ (2050)  , POW (2050)  , SOW  (2050) 

IF  (M .  EQ.  2) GO  TO  20 

PLOT  THE  UNSMOOTHED  POWER  SPECTRUM. 

READ  (7,610)  (ALPHA  (I)  ,1  =  1,20) 

READ  (7,6  11) NF,KA,KB,KC,KA,HB,HC,VA,VB,VC 
READ (7, 612) HT,SYD,SYS 
IF  (KC.GT. 0) GO  TO  10 
WRITE  (1,613)  (FREQ  (I)  ,POW(I)  ,1  =  1, N) 

GO  TO  20 

10  CALL  FLTG  (FREQ,  POW  ,  POW  ,N,  NF ,  KA  ,  KE  ,  KC ,  KC,  H  A,  HB  ,  HC  ,  V  A  ,  VB  ,  VC  ,  AL  PH  A  ,  6  , 
*HT,SYD,SYS,1) 

IF  (M.EQ.  1) GO  TO  40 

PLOT  THE  SMOOTHED  POWER  SPECTRUM. 

20  READ  (7,610)  (ALPHA  (I)  ,1  =  1, 20) 

READ  (7,6 1 1) NF, KA , KB , KC , HA ,HB ,HC , VA,VE,  VC 
READ  (7,612)  HT,SYD, SYS 
IF  (KC. GT. 0) GO  TO  30 
WRITE  (1,6  13)  (FREQ  (I)  , SOW (I)  ,1  =  1 ,N) 

RETURN 

30  CALL  PLTG ( FREQ , SO W , SOW  ,N , NF , KA , KB , KC , KC, HA , HB , HC , V A , VB , VC , ALPHA, 6 , 
*HT, SYD, SYS  ,  1 ) 

40  RETURN 
END 
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Table  B.13*  Subroutine  PLTMEM ,  as  described  in  Chapter  6. 


SUBROUTINE  PLTMEM  (N) 

C  ****:**##****************;( t ************  ********  ********************* 

C  THE  PARAMETERS  TO  BE  READ  IN  THIS  SDBROOTINE  (FROM  UNIT  7)  ARE  THE 

C  SAME  AS  IN  PLTDAT. 

C  IF  NO  PLOT  IS  WANTED  AT  THIS  TIME  BUT  THE  INFORMATION  IS  TO  BE 

C  STORED,  SET  KC  NEGATIVE  AND  ASSIGN  A  FILE  TO  UNIT  2. 

Q  *******************************************  ******  *******  ********** 

620  FORMAT (20 A4) 

621  FORMAT (4I5,6F10.3) 

622  FORMAT (3F10.  3) 

623  FORMAT (8F10. 5) 

*************** ********************************************* ****** 
DIMENSION  ALPHA  (20) 

COMMON  /MEMOUT/  FREQ  (U 097 ) , P (U 09^ ) 

READ  (7,620)  (ALPHA  (I)  , 1  =  1 , 20) 

READ  (7,621)NF,KA,KB,KC,HA,H8,HC,VA,VB,VC 
READ  (7, 6  22) HT,SYD,SYS 
IF  (KC.LE. 0) GO  TO  20 

CALL  PLTG  (FREQ, P,P,N , NF, K A , KB, KC , KC , H A , HB, HC , VA, VB,VC, ALPHA,  6, HT, S 

*  YD, SYS, 1) 

RETURN 

20  WRITE  (2,623)  (FREQ  (I)  , P  (I)  , 1  =  1 , N ) 

RETURN 
END 
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Table  B.m.  Notch  filter  program,  as  described  in  Chapter  6. 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

700 

701 

702 

703 
04 
05 

706 

707 

708 

709 


C 

c 


NOTCH  FILTER 

THIS  PROGRAM  PERFORMS  A  FILTERING  PROCESS  ON  THE  DATA  USING  A 
NOTCH  FILTER  TO  REMOVE  A  PARTICULAR  FREQUENCY. 

IF  REQUIRED,  THE  PROGRAM  WILL  CALCULATE  THE  FREQUENCY  RESPONSE 
FUNCTION;  IT  CAN  ALSO  PLOT  THIS  FUNCTION. 

ZA  -  THE  NEIGHBOURHOOD  OF  THE  UNIT  CIRCLE  DESIRED. 

M  =  NUMBER  OF  TIMES  THAT  THE  FILTER  IS  APPLIED;  ONE  APPLICATION 
INVOLVES  A  FOWARD  AND  A  REVERSE  PASS  OVER  THE  DATA. 

ZB  =  THE  ANGLE  ON  THE  UNIT  CIRCLE  OF  THE  UNWANTED  FREQUENCY,  IN 
RADIANS. 

NFR  =  0,  NO  CALCULATION  OF  THE  FREQUENCY  RESPONSE  FUNCTION. 

=  NUMBER  OF  FREQUENCY  VALUES  REQUIRED  IN  THE  FREQUENCY 
RESPONSE  FUNCTION. 

(POSITIVE)  RESULTS  ARE  PRINTED  AND  PLOTTED. 

(NEGATIVE)  RESULTS  APE  PRINTED  ONLY. 

THE  REMAINING  INPUT  DATA  IS  THE  SAME  AS  FOR  THE  MAIN  PROGRAM  OF 
ANALYSIS. 

THE  PARAMETERS  TO  BE  READ  FOR  PLOTTING  (FROM  UNIT  7)  ARE  THE  SAME 
AS  DESCRIBED  IN  THE  MAIN  ANALYSIS  PROGRAM. 

THIS  PROGRAM  CAN  BE  RUN  IMMEDIATELY  BEFORE  THE  ANALYSIS  PROGRAM  IF 
UNIT  8  IS  ASSIGNED  TO  THE  FILE  FROM  WHICH  THE  ANALYSIS 
PROGRAM  WILL  OBTAIN  ITS  INFORMATION. 

FORMAT  (F 7. 5, 13, FI  0.7, 1  5) 

FORMAT (1315, F5. 1 , F5. 3) 

FORMAT (1H1 , ’THE  COEFFICIENTS  FOR  THE  NOTCH  FILTER  A  RE:  *  , 3 F  1  0 . 7 ) 
FORM  AT  (1H-, 12X, 1  FREQUENCY’ , 1 IX, ' RESPONSE* ,/) 

FO  RM  AT  ( 1  H  ,  1  0  X  ,  F  1  0 . 5 , 1  0  X  ,  F 1  0 . 6 ) 

FORMAT (20 A4) 

FORMAT  (41 5 , 6F1 0 . 3) 

FORMAT  (3  FI  0 . 3) 

FORMAT  (1 5A4) 

FORMAT  ( 1H- , 1  FOR  THE  FREQUENCY  =  ’ ,F10.5,»,  THE  MODULUS  OF  THE  FILT 
*EP  IS  *  , E 1 2 . 5 , *  FOR  EACH  CF  THE', 13,  f  CASCADE  (S)  .  ’ ) 

DIMENSION  X  (41 96)  ,Y  (41 96)  , F  (4  09  7)  , V  (4097)  , ALPHA  (20)  , FORI  (15)  , FOR2  ( 
*15) 

COMPLEX  YA,YB,YC,CMPLX ,CABS 
RF AD  (5,700)  ZA,M,ZB, NFR 
A=1 . 0/ (ZA**2) 

B  =  2 . 0*COS  (ZB) 

C=B/ZA 

READ  (5, 701)  II,  N,  JA,  JB,  JC  ,  JD,  JE,  JF,  JG  ,  JH,  JI,  JJ  ,  JK,  AA,  DT 
WRITE  (6,7  02)  A,  B,  C 
PI2T=2. 0*3. 1 415927*DT 
IF  (NFR.EQ.0)  GO  TO  10 
NU M=I A BS  (NFR) 

M0N=NUM+ 1 

DE  N=  FLOAT (2*NUM) 

WRITE  (6 ,703) 

DO  5  1=1, MUN 
F  (I) =FLO  AT  (1-1) /DEN 


127 


ZC=PI2T*F (I) 

ZD=2 . 0*ZC 

YA=CMPLX  (COS  (ZC)  ,SIN (ZC)  ) 

YB=CMPLX  (COS  (ZD)  ,  SIN  (Z  D)  ) 

YC= (A*  (1 . 0-B*YA+  YB)  ) / ( 1 . 0-C*YA+ A*YB) 

V  (I)  =CA BS  (YC) 

WRITE  (6,704) F(I)  ,  V(I) 

5  CONTINUE 

IF  (NFR. LE. 0)  GO  TO  10 

READ  (7,705)  (ALPHA  (I)  ,1  =  1 , 20) 

READ (7, 706) NF, KA, KB, KC , HA ,HB, HC, VA, VB, VC 
READ  (7, 7  07)  HT, S YD , S YS 

CALL  PLTG (F, V, V , MU N , NF , KA r KB , KC , KC , H A , HB , HC , V A , 7B , VC , ALPH A , 6 , HT, S Y 
*D, SYS, 1) 

10  IF  (II. EQ. 0) GO  TO  15 
READ  (5,708) FORI 
15  READ  (5,7  08) FOR2 
11=11+100 
N0=N+ 1 00 

IF (II. EQ. 100) GO  TO  20 
READ  (5  ,  FOR  1 )  (X  (L)  ,  L=  1 0  1 ,  II) 

20  IM  =  11+ 1 

READ  (5  , FOR2)  (X (L)  , L  =  IM , NO ) 

N 1  =  N+ 1 0  1 
N4=N+ 196 
N6=N+ 1 98 
N7=N+ 199 
N8  =  N  +200 
N9  =  N  +  20  1 
DO  25  1  =  1  , 100 
X(I)  =0.0 
25  Y  (I)  =0.0 

DO  30  I=N 1 , N8 
X  (I)  =0.0 
30  Y  (I) =0.0 
DO  60  L=1,M 
DO  40  1=3, N6 

40  Y  (I)  =A*  (X  (I)  -B*X  (1-1 )  +X  (1-2)  -Y  (1-2)  )  +  C*Y  (1-1) 

DO  50  1=1 , N 4 

50  X(N7-I) =A*  (Y  (N7-I) -B*Y  (N8-I) +Y  (N9-I) -X  (N9-I) )  +  C*X  (N8-I) 

60  CONTINUE 
11  =  0 

WRITE  (8, 7  01)  II,N, JA, JB , JC , JD, JE, JF, JG, JH, JI, JJ, JK, AA,DT 

WRITE  (8, 708)  FOR2 

WRITE (8, FOR2)  (X  (I)  ,1  =  1 01, NO) 

FREQ=ZB/PI2T 

ZD=2.0*ZB 

YA=C  MPLX  (COS  (ZB)  , SIN  (ZB) ) 

YB  =  C MPLX  (COS  (ZD)  ,  SIN  (ZD)  ) 

YC=(A* (1 .0-B*YA+YB) ) /(1.0-C*YA+A*YB) 

YC=YC* YC 
AM  P=CABS (YC) 

WRITE  (6,709)  FREQ, AMP, M 

STOP 

END 
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